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Abstract 

We address the generalized thermodynamics and the collapse of a system of self- 
gravitating Langevin particles exhibiting anomalous diffusion in a space of dimension D. 
This is a basic model of stochastic particles in interaction. The equilibrium states corre- 
spond to polytropic configurations similar to stellar polytropes and polytropic stars. The 
index n of the polytrope is related to the exponent of anomalous diffusion. We consider 
a high-friction limit and reduce the problem to the study of the nonlinear Smoluchowski- 
Poisson system. We show that the associated Lyapunov functional is the Tsallis free 
energy. We discuss in detail the equilibrium phase diagram of self-gravitating polytropes 
as a function of D and n and determine their stability by using turning points arguments 
and analytical methods. When no equilibrium state exists, we investigate self-similar solu- 
tions of the nonlinear Smoluchowski-Poisson system describing the collapse. Our stability 
analysis of polytropic spheres can be used to settle the generalized thermodynamical sta- 
bility of self-gravitating Langevin particles as well as the nonlinear dynamical stability of 
stellar polytropes, polytropic stars and polytropic vortices. 

1 Introduction 

In preceding papers of this series [HI21E1II], we have studied a model of self-gravitating Brownian 
particles enclosed within a spherical box of radius R in a space of dimension D. This model 
can be considered as a prototypical dynamical model of systems with long-range interactions 
possessing a rich thermodynamical structure. For simplicity, we considered the Smoluchowski- 
Poisson (SP) system which is deduced from the Kramers-Poisson (KP) system in a high friction 
limit (or for large times). These Fokker-Planck-Poisson equations were first proposed in pQ as 
a simplified dynamical model of self-gravitating systems. Their relation with thermodynamics 
(first and second principles) was clearly established in terms of a maximum entropy production 
principle (MEPP) and their rich properties (self-organized states or collapse) were described 
qualitatively in this preliminary work. A thorough study of this system of equations was 
undertaken more recently in El HI, complemented by rigorous mathematical results (see 
|3 El HI El M an d references therein) . 

The Smoluchowski equation is a particular Fokker-Planck equation involving a diffusion 
and a drift jlOj . In our model, the drift is directed toward the region of high densities due to 



the gravitational force which is generated by the particles themselves. This retroaction leads 
to a situation of collapse when attraction prevails over diffusion (21 El II]- The SP system also 
provides a simplified model for the chemotactic aggregation of bacterial populations JTj and for 
the formation of large-scale vortices in two-dimensional hydrodynamics ^21 El EJ EH E] • 
It can be shown that the SP system decreases continuously a free energy constructed with the 
Boltzmann entropy |TJ Accordingly, the stationary solutions of the SP system are given 
by the Boltzmann distribution which minimizes Boltzmann's free energy at fixed mass. The 
equilibrium state is then determined by solving the Boltzmann-Poisson (or Emden) equation, 
as in the case of isothermal gaseous stars and isothermal stellar systems (T%1 IT5] . However, 
depending on the value of the temperature T, an equilibrium solution does not always exist and 
the system can undergo a catastrophic collapse. We determined analytically and numerically 
self-similar solutions leading to a finite time singularity (21 El- In OH], we showed that the 
evolution continues after the collapse until a Dirac peak is formed. 

In this paper, we propose to extend our study to a generalized class of Smoluchowski equa- 
tions (201- They can be obtained from the familiar Smoluchowski equation by assuming that 
the diffusion coefficient depends on the density while the drift coefficient is constant (or the 
opposite). They can also be obtained from standard stochastic processes by considering a spe- 
cial form of multiplicative noise [2U 1201 or an extended class of transition probabilities (2*2*1 "2*3*] . 
These equations are consistent with a generalized maximum entropy production principle [20 . 
For simplicity, we shall assume that the diffusion coefficient is a power law of the density. In 
the absence of drift, this would lead to anomalous diffusion. If we take into account a drift 
term and a self-attraction, we have to solve the nonlinear Smoluchowski-Poisson (NSP) system. 
It can be shown [201 121] that the NSP system decreases continuously a free energy associated 
with Tsallis entropy |25j . Accordingly, the stationary solutions of the NSP system are given by 
a polytropic distribution which minimizes Tsallis' free energy at fixed mass. The equilibrium 
state is then determined by solving the Lane- Emden equation, as in the case of polytropic stars 
and stellar polytropes [""HI HH] • Depending on the value of the control parameter and on the 
index n of the polytrope, three situations can occur: (i) the NSP system can relax toward an 
incomplete polytrope maintained by the walls of the confining box (ii) the NSP system can 
relax toward a stable complete polytrope of radius R* < R, unaffected by the box (iii) the NSP 
system can undergo a catastrophic collapse leading to a finite time singularity. 

The paper has two parts that are relatively independent. In the first part of the paper 
(Sees. El and EJ), we study the dynamical stability of stellar systems, gaseous stars and two- 
dimensional (2D) vortices by using a thermodynamical analogy (2111201 123- Stellar polytropes 
maximize Tsallis entropy (considered as a if -function) at fixed mass and energy. This is a 
condition of nonlinear dynamical stability via the Vlasov equation. Polytropic stars minimize 
Tsallis free energy (related to the star energy functional) at fixed mass. This is a condition 
of nonlinear dynamical stability via the Euler- Jeans equations. Polytropic vortices maximize 
Tsallis entropy (considered as a if-function) at fixed circulation and energy. This is a condition 
of nonlinear dynamical stability via the 2D Euler equation. We perform an exhaustive study of 
the structure and stability of polytropic spheres by determining whether they are maxima or 
minima (or saddle points) of Tsallis functional. For sake of generality, we perform our study in 
a space of dimension D. We shall exhibit particular dimensions D = 2, D = 4, D = 2(1 + v^2) 
and D = 10 which play a special role in our problem. The dimension D = 2 is critical because 
the results established for D > 2 cannot be directly extended to D = 2 [3 j. On the other 
hand, the nature of the caloric curve changes for D = 4 and D = 10. This extends the study 
performed by Taruya & Sakagami (23 12B] and Chavanis (2H1 12EJ for D = 3. 

In the second part of the paper (Sec. EJ), we study the dynamics and thermodynamics of self- 
gravitating Langevin particles experiencing anomalous diffusion. Their equilibrium distribution 
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minimizes Tsallis free energy at fixed mass. This is a condition of thermodynamical stability 
in a generalized sense. This is also a condition of linear dynamical stability via generalized 
Fokker-Planck equations (in the present context the NSP system) fZ\ I24[ l2*Uj . Thus, the stability 
analysis of Sec. El can also be used in that context. When the static solution does not exist or is 
unstable, the system undergoes a catastrophic collapse. In Sec. HI we show that the NSP system 
admits self-similar solutions describing the collapse and leading to a finite time singularity. The 
density decreases at large distances as p ~ r~ a . In the canonical situation (fixed T), the scaling 
exponent is a n = 2n/(n — 1) where n is the polytropic index. We also consider a microcanonical 
situation in which the generalized temperature T(t) varies in time so as to rigorously conserve 
energy. In that case, the scaling equation has solutions for a n < a < a max (n,D) where 
a max (n, D) is a non-trivial exponent. The value of a effectively selected by the system is 
determined by the dynamics. In Sec. we perform direct numerical simulations of the SP and 
NSP systems with higher accuracy than in Ref . [2] . We confirm the scaling regime and discuss 
the value of the scaling exponent. In the microcanonical situation, we give numerical evidence 
that a = a n so that the temperature T(t) is finite at the collapse time (this observation has 
been made independently by [9 for the SP system). However, the convergence to this value 
is so slow that the evolution displays a pseudo-scaling regime with a n < a < a max for the 
densities achieved. We explain the physical reason for this behavior and we conjecture that 
a ~ ct max will be reached in more realistic models with non-uniform temperature [20] • This 
paper closely follows the style and presentation of our companion paper for isothermal spheres 
P|. These two papers complete the classical monograph of Chandrasekhar on self-gravitating 
isothermal and polytropic spheres in D = 3 [TSj . 



2 Dynamical stability of systems with long-range inter- 
actions 

2.1 Stellar systems 

Let us consider a collection of N stars with mass m in gravitational interaction. They form a 
Hamiltonian iV-body system with long-range (Newtonian) interactions. We work in a space of 
dimension D and enclose the system within a spherical box of radius R. Let f(r,v,t) denote 
the distribution function of the system, i.e. f(r,v,t)d D rd D v gives the mass of stars whose 
position and velocity are in the cell (r, v; r + d D r, v + d D v) at time t. The integral of / over 
the velocity determines the spatial density 



(1) P = Jfd D v. 
The total mass of the configuration is 

(2) M = J pd D r. 

In the mean-field approximation, the total energy of the system can be expressed as 

(3) E = \j f ^ dDrd ° V + \ J P$ d ° r = K + W > 

where K is the kinetic energy and W the potential energy. The gravitational potential $ is 
related to the density by the Newton-Poisson equation 

(4) A$ = S D Gp, 
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where So is the surface of a unit sphere in a D-dimensional space and G is the constant of 
gravity. 

For fixed N ^ 1 and t — > +00, the system is expected to reach a statistical equilibrium state 
described by the classical Boltzmann entropy Sb [/] = —/ / In fd D rd D v. However, the relax- 
ation time t re i ax due to "collisions" (more properly close encounters) is in general considerably 
larger than the dynamical time to so that this statistical equilibrium state is often not physi- 
cally relevant [HUES- This is the case in particular for elliptical galaxies where t re i ax ~ ^777/0 
with iV ~ 10 12 , while their age is ~ lOOt^ [TH]. For t <C t re i ax and N — > +00, the dynamics of 
stars is described by the Vlasov equation 

where F = — V$ is the gravitational force determined by the Poisson equation (@J). For any 
if-function 



i D rd D 



(6) S[f] = - J C{f)d l 

where C is convex, i.e. C" > 0, it can be shown that the variational problem 

(7) Max S[f] at fixed E[f],M[f], 

determines a stationary solution of the Vlasov equation with strong (nonlinear) dynamical 
stability properties [^SlIHlj. Such solutions can result from a process of (possibly incomplete) 
violent relaxation [2E] ■ Introducing Lagrange multipliers, the first order variations 5S — (35 E — 
a6M = lead to 



(8) &{f) = -?[- + * 



v 2 



a. 



Therefore, / = /(e) where e = y + $ is the energy of a star by unit of mass. These distribution 
functions, depending only on the energy, form a particular class of stationary solutions of the 
Vlasov equation. Other solutions can be constructed with the Jeans theorem JH] but their 
stability is more difficult to investigate. The conservation of angular momentum can be easily 
included in the foregoing discussion 



2.2 Barotropic stars 

Let us now consider a self-gravitating gaseous system described by the Euler- Jeans equations 
(9) ^ + V( P u) = 0, 



(10) — + (u-V)u = — Vp-V$. 

at p 

We assume that the gas is barotropic with an equation of state p = p(p). The most important 
examples of barotropic fluids are those that are isentropic or adiabatic, that is those whose 
specific entropy is constant. If s = Cst., the first principle of thermodynamics du = —pdv + Tds 
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(where v = 1/p) reduces to du = j^dp, where u is the internal energy by unit of mass. The 
total energy of the fluid is therefore 

(11) W[p] = J pJ^dp'd D r+ 1 - f p^d D r + f p^d D r. 

The first term is the internal energy, the second the gravitational energy and the third the 
kinetic energy associated with the mean motion. As argued in the variational problem 

(12) Min W[p] affixed M\p], 

determines a stationary solution of the Euler- Jeans equations with strong (nonlinear) dynamical 
stability properties. We shall not prove this result here but it is expected to follow from 
relatively standard methods of stability theory. The solutions of this variational problem satisfy 
the condition of hydrostatic balance 

(13) Vp=-pV$, 

between pressure and gravity. The foregoing results can be extended to barotropic stars rotating 
rigidly with angular velocity 12 [35] . 



2.3 Two-dimensional vortices 

Let us finally consider a collection of N point vortices with circulation 7 in 2D hydrodynamics. 
They form a Hamiltonian system with long-range (logarithmic) interactions. We call u the 
vorticity, ip the stream function and u = —zxVip the velocity field. We also note r = j ud 2 r 
the circulation and E = | J uipd 2 r the energy. For fixed N ^ 1 and t — > +00, the system is 
expected to reach a statistical equilibrium state described by the classical Boltzmann entropy 
Sb[u] = — f uj\n.ud 2 Y [36J. However, the relaxation time t re i ax due to "collisions" is in general 
considerably larger than the dynamical time to so that this statistical equilibrium state is often 
not physically relevant [IH1 EE]. For t t Te i ax and iV — > +00, the dynamics of point vortices is 
described by the 2D Euler-Poisson system 

(14) ^ + u . V cu = 0, 



(15) uj = -Alp. 

The Euler equation also governs the dynamics of 2D incompressible and inviscid continuous 
vorticity fields. For any if-function 



(16) S[co] = - J C(uj)d 2 r } 

where C is convex, i.e. C" > 0, it can be shown that the variational problem 

(17) Max S[u] affixed E[u],T[u], 

determines a stationary solution of the 2D Euler equation with strong (nonlinear) dynamical 
stability properties jHZj- Such solutions can result from a process of (possibly incomplete) 
violent relaxation [TJj. Introducing Lagrange multipliers, the condition 5S — (35 E — a5T = 
leads to 

(18) c'H = - «■ 

Therefore, u = uj(ip), which is the general form of stationary solutions of the 2D Euler equation 
for domains with no specific symmetries. The conservation of angular momentum and impulse 
can be easily included in the foregoing discussion [20j. 
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2.4 Thermodynamical analogy 

Since the variational problems (J7J), (fT2*j) and (|17|) are similar to the usual variational problems 
that arise in thermodynamics (with the Boltzmann entropy), we can develop a thermodynamical 
analogy to analyze the dynamical stability of stellar systems, gaseous stars and 2D vortices 
j2Hl 1201 1231 HZ! • I n this analogy, the functional S plays the role of a generalized entropy, (3 is a 
generalized inverse temperature, (3(E) a generalized caloric curve etc... The variational problem 
(j7J) is similar to a condition of microcanonical stability. We can also introduce a generalized 
free energy F[f] = E[f] — TS[f] which is the Legendre transform of S[f}. The minimization of 
F[f] at fixed T and M[f] is similar to a condition of canonical stability. This is equivalent to 
first minimize F[f] at fixed p(r) to get /*(r, v) and then to minimize F[p] = F[f*], calculated 
with /*, at fixed M[p\. Now, it can be shown j2Hj that F[p] is precisely the functional (JTTJ) 
with u = 0. Therefore, the variational problem (|12|) is similar to a condition of canonical 
stability. Since canonical stability implies microcanonical stability (but not the converse) [20 , 
we conclude that "stellar systems are stable whenever corresponding barotropic stars are stable" 
which provides a new interpretation of Antonov's first law [2*o] . 

In 2D hydrodynamics, the variational problem (|17|) is similar to a condition of microcanon- 
ical stability. It is stronger than the maximization of J[u] = S[u] — (3E[u] at fixed (3 and T[u] 
(canonical stability), which is just a sufficient condition of nonlinear dynamical stability. It 
is not a necessary condition of stability if the ensembles are inequivalent (i.e. if the "caloric 
curve" presents bifurcations or turning points). The Arnold's theorems just provide sufficient 
conditions of canonical stability (see, e.g., |2U|). Therefore, in the domain of inequivalence (cor- 
responding to a region of "negative specific heats" ) a flow can be nonlinearly dynamically stable 
while it violates Arnold's theorems. This has important implications in jovian fluid dynamics 

maun. 

The preceding arguments also apply to other systems with long-range interactions such as 
the HMF model for example • This opens the route to many generalizations by changing 
the potential of interaction and the "generalized entropy" (H-function). In this paper, we shall 
specialize in the case of particles interacting via a Newtonian potential (e.g., self-gravitating 
systems, 2D vortices,...). We shall also consider a special form of if-function, known as Tsallis 
entropy, leading to power-laws distributions (polytropes). 

3 Equilibrium structure of polytropic spheres in dimen- 
sion D 

3.1 Stellar polytropes 

Let us consider a particular class of stationary solutions of the Vlasov equation called stellar 
polytropes [19] . Nonlinearly dynamically stable solutions maximize the Tsallis entropy 

(19) Sq = ~^l J(f q ~f) dDrdD V' 

at fixed mass M and energy E, where q is a real number. For q — > 1, Eq. ([TH]) reduces to 
the ordinary Boltzmann entropy describing isothermal stellar systems. We emphasize that, 
in the present context, Tsallis and Boltzmann entropies are particular iJ-functions (not true 
entropies) that are related to the dynamics, not to the thermodynamics. Still, due to the 
thermodynamical analogy discussed in Sec. 12.41 we shall use a thermodynamical langage to 
study the dynamical stability problem. In this analogy, the dynamical stability criterion for 
stellar polytropes corresponds to a microcanonical stability condition. 
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The critical points of entropy at fixed mass and energy satisfy the condition 



(20) 



5S q - (35E - X5M = 0, 



where (3 = 1/T and A are Lagrange multipliers (T is the temperature and A the chemical poten- 
tial in the generalized sense). The variational principle (j2*Uj) leads to the polytropic distribution 
function 



(21) 



/(r,v) 



(9-1)0 
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— + $(r) 



9-1 



where /i = [1 — (g — l)X]/q. We define the polytropic index n by the relation 

D 1 
~2 + g- 1" 



(22) 



n 



We first consider the case (g — l)0/g > and allow to take negative values. Then, Eq. 
(12 lj) can be rewritten 



(23) 



where 



(24) 



.4 



9 



a-d>- y 



9-1 



l-(g-l)A 
(9-1)0 ' 



If n > D/2 (i.e., q > 1 and > 0), /'(e) < where e = ^- + $ is the stellar energy. 
Therefore, high energy particles are less probable than low energy particles, which corresponds 
to the physical situation. Equation (|2*3j) is valid for v < v max = \/l{a — Wj. If v > v max = 
y2(a — $), we set / = 0. This distribution function describes stellar polytropes which were 
first introduced by Plummer If n = D/2 (i.e., g — > oo), the distribution /(e) is a step 

function. This corresponds to the self-gravitating Fermi gas at zero temperature (white dwarfs). 
In D = 3, classical white dwarf stars are equivalent to polytropes with index n = 3/2 jTHJ. In 
/^-dimensions, classical "white dwarf stars" are equivalent to polytropes with index 



(25) 



^3/2 



D 
~2 ' 



If n — > +oo (i.e. q — > 1), we recover isothermal distibution functions. If n < D/2 (i.e., g < 1 
and g/3 < 0), high energy particles are more probable than low energy particles: /'(e) > 0. This 
situation is unphysical but it can be considered at a formal level. The distribution function 
diverges like (1 — v /v max ) n ~ D l 2 as v — > v max . The moments of / converge if and only if 
n > D/2 — 1. Therefore, if D/2 — 1 < n < D/2 (i.e., q < and > 0), stellar polytropes exist 
mathematically but they are not physical. If n < D/2 — 1 (i.e., < q < 1 and < 0), stellar 
polytropes do not exist. 

We now consider the case (g — l)0/g < 0. Then, Eq. (|21|) can be rewritten 



(26) 



f = A 



a + $ + 



2 



where 
(27) 



.4 



'(1-9)0" 



9-1 



a: 



l-(g-l)A 
(1-9)0 ' 
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If n > D/2 (i.e., q > 1 and (3 < 0), the model is ill-posed because f(v) diverges for v — > +oo. 
If n < D/2, the distribution function goes to zero like v-^ D - 2n ^ for v +oo. If < n < D/2 
(i.e., q < 1 - 2/D), the density p = / fd D v does not exist. If -1 < n < (i.e., 1 - 2/D < 
q < D/(D + 2)), the density exists but not the pressure p = jj f fv 2 d D v. If n < — 1 (i.e., 
D/(D + 2) < q < 1 and (3 > 0), the density and the pressure exist. 

In conclusion, only positive temperatures states are physical. If n > D/2 (i.e. q > 1), the 
system is described by the distribution function (}2l?j) . If n < — 1 (i.e., D/[D + 2) < g < 1), 
the system is described by the distribution function ()26|) . In this paper, we shall only consider 
stellar polytropes with index n > D/2. Using Eq. (j2HJ), the spatial density p = J fd D v and 
the pressure P = j) f fv 2 d D w can be expressed as 

(28) p = 2 D/2 - 1 AS D {a - ®) n B(D/2, n + 1- D/2), 



(29) p = ^^2 D/2 - 1 AS D (a - ^) n+1 B(D/2, n + 1- D/2), 

with B(a, b) being the beta-function. In obtaining Eq. (j2HJ), we have used the identity B(m + 
1, n) — mB(m, n)/(m + n). Eliminating the gravitational potential between these two relations, 
we recover the well-known fact that stellar polytropes satisfy the equation of state 

(30) p = Kp\ 7 = 1 + ^, 

like gaseous polytropes (see below). In the present context, the polytropic constant is given by 

(31) K = ——h D ^- 1 S D AB(D/2,n+ 1 - D/2) j . 



Using Eqs. (jzlj) and (|28j) . the distribution function (j25j) can be written as a function of the 
density as 

v 2 /2 i"- D/2 



(32) / ' 



p l/n 



(n + l)K 



Z 
with 

(33) Z = 2 D/2 - 1 S D B(D/2, n + 1- D/2) [K(n + 1)} D/2 . 

The distribution function (j32j) can also be obtained by maximizing S q [f] at fixed M, E and 
p(r), or equivalently at fixed K = | j fv 2 d D rd D \ and p(r). It is then possible to express the 
energy ® and the entropy (fT9"j) in terms of p(r) and T. Using Eq. (J3~2~j) . it is easy to show that 

(34) E = ^J pdDr+ \ j P® d ° T > 

(35) S q [p] = - L - J P d D r - M 

In arriving at Eq. ()35|). we have used the identity Bim, n + 1) = nB(m,n)/(m + n). Proceeding 
carefully, we can check that for q — > 1, Eq. (j33j) reduces to the Boltzmann entropy expressed in 
terms of hydrodynamical variables (see Eq. (9) of jS])- The problem of the stability of stellar 
polytropes now amounts to determining maxima of S q [p] at fixed E[p] and M[p]. 
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3.2 Gaseous polytropes 



We shall consider a particular class of barotropic stars called gaseous polytropes. They are 
characterized by an equation of state of the form 

(36) p = Kp\ 7 = 1 + -, 

n 

where K is a constant. We recall that gaseous polytropes are described by a local thermody- 
namical equilibrium condition and not by a distribution function of the form (j21j) . except in 
the isothermal case n — > +oo (see [29 J. Their energy is 

(37) W[p) = n j pd D r p$d D r + / y^r. 

Using Eqs. (J3~2|) and (J33j) . we note that the free energy of stellar polytropes F q = E — TS q is 

(38) F q [p] = l -J p$d D v + n J pd D r, 

within an additional constant. Taking u = 0, we check on this specific example that F q [p] = 
W[p]. In fact, this relation is general as shown in Therefore, the dynamical stability 

of gaseous polytropes can be settled by studying the minimization of F q [p] at fixed mass. In 
the thermodynamical analogy of Sec. 121 this corresponds to a canonical description. We note 
finally that the free energy ()38|) can be written explicitly 



(39) F q [p] = l - J p$d D r + J (P 7 ~~ P)d D r. 

This can be viewed as a free energy F 7 = E — KS 1 associated with a Tsallis entropy in position 
space S 1 = — / (p 7 — p)d D r, where 7 plays the role of the q parameter and K the role of a 
temperature. The parameters q and 7 are related to each other by 7 = l + 2(g — l)/[2+D(q — 1)]. 
The equilibrium distribution can be written 



(40) p 



7-1 



7-1 

which is equivalent to Eq. (j2HI)- When considering gaseous polytropes, we shall allow for 
arbitrary value of the index n > 0. 



3.3 The D-dimensional Lane-Emden equation 

The configuration of a stellar polytrope is obtained by substituting Eq. (J2l?j) in the Poisson 
equation (HJ) using Eq. (JU). This yields a self-consistent mean-field equation for the gravitational 
potential $. An equivalent equation can be obtained by substituting the equation of state (|3T?|) 
in the condition of hydrostatic equilibrium (fE^) . as for gaseous polytropes (the equivalence 
between these two approaches is shown in |2E])- Using the Gauss theorem 

rf$ GM(r) 



D-l 



dr r 

where M(r) = Jq pSrjr' D ~ 1 dr' is the mass within the sphere of radius r, we can rewrite Eq. 
(fT3*j) in the form 

tAn\ 1 d ( rD ~ ld P\ c n 

r u dr \ p dr J 
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which is the fundamental equation of hydrostatic equilibrium in D-dimensions. For the poly- 
tropic equation of state (|5U|). we have 



(43) 



K(n+1)- 



1 d ( D _ 1 dp 



.d-i dr 



l/n 



dr 



-SoGp. 



The case of isothermal spheres with an equation of state p = pT is recovered in the limit 
n — ► +00. To determine the structure of polytropic spheres, we set 



(44) 



P = Po9\ £ 



SdGpq 



1-1/n-i 1/2 



K(n + 1) 

where po is the central density. Then, Eq. (|4"5|) can be put in the form 



(45) 



1 d ( dD _ x d6 y 



which is the D-dimensional generalization of the Lane-Emden equation [18 j . For D > 2 and 

D 



(46) 



n > 



D-2 



n 3 , 



Eq. (|43j) has a simple explicit solution, the singular polytropic sphere 

(47) e s -- 



2[{D - 2)n - D}) ~ 

£ n-l. 



(n- l) 2 

The regular solutions of Eq. (|4~HJ) satisfying the boundary conditions 
(48) 9 = 1, 9' = at f = 0, 

must be computed numerically. For £ — > 0, we can expand the solutions in Taylor series and 
we find that 



(49) 



9=1 — 777^^ 



n 



2/T 8D(D + 2) 



e 4 + ... 



To obtain the asymptotic behavior of the solutions for £ — > +00, we note that the transformation 
t = ln£, 9 = £- 2 /(™-i) z changes Eq. (jSJ) in 



(50) 



d 2 z (£> - 2)n - (D + 2) 
cit 2 n — 1 (it 



2[L>+ (2-D)n] 
(n-l) 2 



For £> < 2 or for D > 2 and 
(51) 



D + 2 
n< D-2= n5 > 



the density falls off to zero at a finite radius R*. This defines a complete polytrope of radius 
i?*. If we denote by £1 the value of the normalized distance at which 9 = 0, then, for £ — > £1, 
we have 



(52) 



= -&0i 



+ 



+ 



+ ... 
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On the other hand, for D > 2 and n > n 5 , Eq. (J5(J|) corresponds to the damped motion of a 
fictitious particle in a potential 

/ x D + (2-D)n 9 1 „,i 

(53) V(z) = 1 J z 2 + -— z n+1 , 

[n — iy n+1 

where z plays the role of position and t the role of time. For t — > +oo, the particle will come 
at rest at the bottom of the well at position z = {2[(D — 2)n — D]/{n — l) 2 } 1 ^ 1 . Returning 
to original variables, we find that 

(54) ^{*^!}V=»„ for ^ + oo. 

Therefore, the regular solutions of the Lane-Emden equation ()45|) behave like the singular 
solution for £ — > +oo. To determine the next order correction, we set z = z + z' with z' <C 1. 
Keeping only terms that are linear in z', Eq. (J5U|) becomes 

d 2 z> (D - 2)n - (D + 2) cfe' 2[(D - 2)n - D] , _ 
(55j rft^ + n^l dt + W=l " 

The discriminant of the second order polynomial associated with this equation is 

_ —(D - 2)(10 - P)n 2 - 2(D 2 - 8D + 4)n + QD - 2) 2 
(5bj A(nj - ^— jy 2 . 

For n — > +oo, A(n) ~ — (D — 2) (10 — -D). Furthermore, A(n) = for 



-D 2 + 8L> - 4 ± 8VD^T 
(57) n± = (D-2)(10-D) • 

For 2 < £) < 10, it is straightforward to check that n_ < n + < 725. Therefore, for n > n^, A 
has the sign of —(D — 2) (10 — .D) which is negative. Thus 



(58) = s ji + ^ cos ^^l n £ + ^j, (e^+00), 

where 6 = [(D — 2)n — (D + 2)]/(n — 1). The density profile intersects the singular 
solution (}4Tj) infinitely often at points that asymptotically increase geometrically in the ratio 
1 : exp{27r/v / — A} (see, e.g., Fig. 1 of Ref. [2E] for D = 3). For D > 10, we have n + < n 5 < n_. 
Therefore, if < n < n_, A has the sign of (D — 2) (10 — D) which is negative and the 
asymptotic behavior of the solutions is still given by Eq. ()58|). However, for n > n_, A(n) is 
positive and therefore 

(59) o = oIi + 1(a£# + BZ=P)\, (C^+oo). 



e /2 

Finally, for n = A = and we have 

(60) # = ^i + _L(Aln£ + 5)}, (e^+00). 

For D = 10, n _ ^ +00 so that the asymptotic behavior of 9 is given by Eq. (|5K|) for n < +00. 
Since 9™ ~ £-2n/(n-i) a ^ } ar g e distances, the configurations described previously have an "infinite 
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mass" which is clearly unphysical. In the following, we shall confine these configurations to a 
"box" of radius R as for isothermal spheres [3]. Such configurations will be called incomplete 
polytropes. We note that the self-gravitating Fermi gas at zero temperature forms a complete 
polytrope only if n 3 / 2 < n 5 , i.e. D < 2(1 + y/2) ~ 4.83. 

The Lane-Emden equation can be solved analytically for some particular values of the 
polytropic index. For n = 0, which corresponds to a body with constant density, we have 

(61) 6 = 1- 6 = V2D. 

For n = 1, Eq. (J4*5*|) reduces to the D-dimensional Helmholtz equation. For D = 1, 

(62) fl = cos£, £i = f- 
For D = 2, 

(63) = •/<>(£), 6= 2.40482... 
Finally, for D = 3, performing the change of variables = x/£, we get 

(64) 6 .,. 

The Lane-Emden equation can also be solved analytically in any dimension of space D > 2 for 
the particular index value n 5 . The solution is 

(65) 9 5 ' 



^ X ^ D(D—2) ) 



D-2 1 

D(D-2) , 

as can be checked by a direct substitution in Eq. (JJSJ). For = 3, we recover the Schuster 
solution ^H]. We note that 9$ ~ £ 2 - D for £ — > +oo implying a finite mass. This contrasts with 
the asymptotic behavior (|54|) of the solutions of the Lane-Emden equation with index n > 775. 
For D = 2, ri5 ~~ * +oo and we are lead back to the isothermal case where an analytical solution 
is also known 

Finally, for D = 1, the Lane-Emden equation reduces to the form 

This equation corresponds to the motion of a fictitious particle in a potential V(8) = 9 n+1 / (n + 
1), where 9 plays the role of position and £ the role of time. The first integral of motion is 



(66) 



(67) £ = I«V" +1 



2\d£j 77 + 1 

The "energy" E is determined by the boundary condition pKj) yielding E = l/(n + 1). Thus, 
the solution #(£) can be written in integral form as 

f 1 dx ( 2 \ 1/2 t 

(68) X (i-^/ 2 -UnJ e - 

Except for n = and n = 1, it seems not possible to obtain in a closed form. However, 
the zero £i of is explicitly given by 



(69) 6 



n + i\ 1/2 ^r(gs 



7T- 



' 3+ra ' 
>2(l+n)< 



Therefore, £i ~ \fnj2 for 77 — » +00. Furthermore, according to Eq. (|67|1. we have #'(£i) 
-(2/(r7 + l)) 1 /2. 
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3.4 The Milne variables 



As is well-known jTH], polytropic spheres satisfy a homology theorem: if is a solution of the 
Lane-Emden equation, then A 2 ^ n ~^9(A^) is also a solution, with A an arbitrary constant. This 
means that the profile of a polytropic configuration of index n is always the same (characterized 
intrinsically by the function 9), provided that the central density and the typical radius are 
rescaled appropriately. Because of this homology theorem, the second order differential equation 
(|4"5Jl can be reduced to a first order differential equation for the Milne variables 

(70) u = — — and v = — — . 

Taking the logarithmic derivative of u and v with respect to £ and using Eq. (|45jl . we get 

1 du 1 , . 
71 - — = -{D-nv-u , 

udt, 4 



(72) -^ = -(2-D + u + v). 

vd£ £ 

Taking the ratio of the foregoing equations, we obtain 

udv u + v — D + 2 



(73) 



v du u + nv — D 



The solution curve in the (u, v) plane is plotted in Figs. ^IH for different values of D and n. 
The (u,v) curve is parameterized by £. It starts, at £ = 0, from the point (u,v) = (D,0) 
with a slope (dv/du)o = —{D + 2)/nD. The points of horizontal tangent are determined by 
u + v — D + 2 = and the points of vertical tangent by u + nv — D = 0. These two lines 
intersect at 

(74) u, = i D - 2 >- D - 2 



n — 1 n — 1 

which corresponds to the singular solution (|47|). 

The Milne variables can be expressed in terms of £ explicitly for particular values of the 
polytropic index. For n — 0, using Eq. (|61|h we have 



(75) u — D, 
For n = n 5 , using Eq. (j65|) . we have 

(76) tx D 



2£ 5 



1 + ^1 + 



£ 2 



D(D-2) x 1 D(D-2) 

Eliminating £ between these two relations, we get 

(77) -^v + u = D. 

More generally, using the asymptotic behavior of determined in Sec. 13.31 we can 

deduce the form of the solution curve in the (u,v) plane. For (2 < D < 10, n > n 5 ) and for 
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Figure 1: Phase portrait of the Lane-Emden equation in the (u, v) plane for D < 2 (specifically 
D = 1). The value of the polytropic index is indicated on each curve. For n — > +oo, denoted 
^isoth, we recover the phase portrait of isothermal spheres. 



{D > 10, n 5 < n < n_), the solution curve spirals indefinitely around the limit point (u s ,v s ). 
For (D > 10, n > the curve reaches the point (u s ,v s ) without spiraling. For D = 10, 
n_ — > +oo and n 5 = 3/2. For D < 2 and for (D > 2, n < 715), the (u, v) curve is monotonic 
and tends to (u, v) = (0, +00) as £ — > £1. More precisely, using Eq. we have for n < +00 

n+l 

For n — > +00, £1 — > +00 and we are lead to our previous study jS]. For D = 2, n 5 — » +00. 

For n — > +00, it is easy to check that the Lane-Emden function (for polytropes) is related 
to the Emden function ip (for isothermal spheres) by the equivalent 

(79) 6 n ~ -e"^. 

n 

This suggests to introduce the variables U = u and V = (n + l)v instead of (J7UJ). For n — > +00, 
U and V tend to the Milne variables u, v defined in the case of isothermal spheres (see [Sj). We 
could have introduced these variables since the beginning but we prefer to respect the notations 
used by Chandrasekhar in his classical monograph [TS] . 



3.5 The thermodynamical parameters 

For an incomplete polytrope confined within a box of radius R, the solution of Eq. (|45|) is 
terminated by the box at the normalized radius (see Eq. (|4*4*|l ) 



(80) 



a 



S D Gpl~ l/n 
K(n + 1) 



1/2 



R. 
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Figure 2: Phase portrait of the Lane-Emden equation in the (u, v) plane for D = 2. At 
this dimension n 5 — > +00, so that the 2D-Schuster solution (J55)l becomes equivalent to the 
2.D-isothermal solution PJ. In the (u, v) plane this corresponds to a straight line. 




Figure 3: Phase portrait of the Lane-Emden equation in the (u,v) plane for 2 < D < 10 
(specifically D = 3). 
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Figure 4: Phase portrait of the Lane-Emden equation in the (u, v) plane for D > 10 (specifically 
D = 15). For D = 10, n_ — > +oo. Therefore, for n > = 3/2, the phase portrait is always a 
spiral except for the index n = +oo for which the spiral is reduced to a point. 




We shall now relate the parameter a to the polytropic constant K (or generalized temperature 
T) and to the energy E. Starting from the relation 



(81) 



M 



R 



pS D r D 1 dr = S D p 



K(l + n) ] D/2 f a 



,SdGp 

and using the Lane-Emden equation we get 



l-l/n 



(82) 



M 



-SdPo 



K(l + n) 



l-l/n 



D/2 



a D - l 6'(a). 



Expressing the central density in terms of a, using Eq. ()80|) . we obtain after some rearrange- 
ments 



(83) 



M 



-S 



D 



K(l + n) 



n-l (D-2)n-D n + 1 

i? „_i a* 1 - 1 0(a). 



For a complete polytrope with radius < R, we need to stop the integration at £ = £i. Thus, 
the equivalent of the foregoing equation is the "mass-radius" relation 



(84) 



(Q-2)(n 3 -») im+Ti) »=i 

M~R„ n = v — , 



l/n 



where o; ra is defined by Eq. f!78|) . For n = n^, the mass is independent on the radius. This 
mathematical property is related to the limiting mass of Chandrasekhar for relativistic white 
dwarf stars ^Bj- For n = 1, the radius is independent on mass. For incomplete polytropes, the 
parameter 



(85) 



_ M 
Sd 



S jjG 



K(l + n) 



1 



(Q-2)n-Q ' 

R ™-i 
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can be considered as a normalized inverse temperature [29] . Indeed, for a given mass M and 
box radius R, it simply depends on the polytropic constant K which is itself related to (3 via 
Eqs. (|31|) and (|24jl . In addition, for n — > +00, the parameter 77 reduces to the corresponding 

(3GMm 

R D-2 , 

Eq. (|55j) can be rewritten 



one for isothermal spheres (77 ~ rj^/n, 77^ = /3 = 1/kT) In terms of this parameter, 



(86) 77 = -a£r0'(a). 

This relation can be expressed in terms of the values of the Milne variables at the normalized 
box radius. Writing uq = u(a>) and vq = v(a), we get 



(87) 77 = («o< 



n-l 



For a given box radius R and a given polytropic constant (or generalized temperature T), this 
equation determines the relation between the mass M and the central density po (through the 
parameter a). For D < 2 and for (£> > 2, n < 77.5), the normalized box radius a is necessarily 
restricted by the inequality a < £1. For the limiting value a = £1, corresponding to a complete 
polytrope with radius i?* = i?, we have 

(88) 77(6) = w n , 

More generally, for complete polytropes with radius i?* < R, we have 



(89) 77 = u r , 



n(D-2)-D 

R* 
~R 



Coming back to incomplete polytropes, we note that for the index n = n 5 , Eq. (|H7|) can be 
written explicitly as 



D+2 

a 2 

(90) r] 



For a — > +00, we observe that 77 ~ (3/a) 1 / 2 for D = 3. 

The computation of the energy is a little more intricate. We first recall the expression of 
the Virial theorem in dimension (see Appendix lB|) : 

(91) 2K + (D -2)W = DV D R D p(R), 

where Vd = Sd/D is the volume of a hypersphere with unit radius. For D = 2, the expression 
of the Virial theorem is (see Appendix lB|) : 

(92) 2AT- = 2irR 2 p(R). 

The potential energy of a polytrope can be calculated as follows. Combining the condition of 
hydrostatic equilibrium (|T3*j) with the equation of state (jHUJ), we get 

This equation integrates to give 

(94) , = _(„ + 1 )g-g|) +.(*). 
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Inserting this relation in the integral (j2J) defining the potential energy W and recalling that the 
kinetic energy can be written K = (D/2) J pd D r, we obtain 

(95) W = -I(n + l)K + J (n + 1)^||m + ^M$(,R). 

For D ^ 2, $(#) = -GM/[(D -2)R°- 2 ] and for L> = 2, we take the convention $(R) = (see 
Appendix EJ. Eliminating the kinetic energy between Eqs. JSJ), (jHSJ) and (JHHJl . we obtain for 

(96) ^ = 7^ ^ -|(n+l)V D i? D j9( J R)-(n + l)4^M+ G '' U " 



D + 2- (D-2)n i v ' ^ y ' v p(i?) (D-2)R D ~ 
and for D = 2: 

(97) W = -(n + 1)^ - \(n + l)nR 2 p(R) + \(n + 1)^^. 

For complete polytropes for which p(R*)/ p(R*) = 0, we obtain the D-dimensional generalization 
of the Betti-Ritter formula 



—D GM 2 
D + 2 — (D — 2)n (D - 2)R? 



( 98 ) W = —n , n ~^-2 ( D * 2 )' 



(99) 



W 



-{n + l)^f + \GM^J^) (D = 2). 



We note that for D > 2, the potential energy is infinite for n = n 5 (while the mass is finite). 
Returning to incomplete polytropes, the total energy E = K + W can be written for D^2: 



E 



-1 



D(A-D) 
D + 2 - (£> - 2)n \ 2(D - 2) 



GM 2 



(n + l)(D-2)M 



Kg) 
p(i2) 



(100) 

and for D 
(101) 



-DV D R D (n + 1 - D)p(fl) 



2: 



E = -(n- 1) 



GM 2 



l -(n-l)^p(R) + \(n+l) P j^M. 



Expressing the pressure p(R) in terms of the Lane-Emden function 9(a) using Eqs. fl^OI) and 
flSJ), using Eqs. (fKUJ) and (jHHJ) to eliminate the central density po an d the polytropic constant 
K (or temperature T), and introducing the Milne variables (|70|). we finally obtain for D ^ 2: 



(102) A = - 

and for D = 2: 
(103) 



£>-2 



GM 2 



(D - 2)n -(D + 2) 



D(A — D) 
_2(D-2) 



D-2\ n+l-Duo 

V Q J 71 + 1 W 



A = -in — 1 J H 

8 v ; 4 n + 1 v 



1 
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For D < 2 and for (D > 2, n < n 5 ), the normalized box radius a in necessarily restricted by 
the inequality a < £1. For the limiting value a = £1, corresponding to a complete polytrope 
with radius R* = i£, we have 

(104) A(£i) = A„, 
with 

(105) A„ = -£>(4-D) 

v ; 2(D-2)[(D-2)n-(D + 2)] V ; 



(106) A n = ^(ra-1) (D = 2). 

o 

More generally, for complete polytropes with radius < i?, the dimensionless energy is 

R ,D-2 



(107) A = A n ( — ) , (1^2), 

ft* 



(108) A = ±(n-l)-±ln(:|) (D = 2) - 
Eliminating i2* between Eqs. flH^j) . ()107|) and (|108|) . we obtain 

(n-l)(D-2) (n-l)(P-2) 

(109) At? = \ n Un {D ~ 2) ~ D , (-D ^ 2), 



(HO) 



A = -in - I) 



2 In 



( V_ 



(D = 2). 



This defines the branch of complete polytropes in the A — rj plane. Coming back to incomplete 
polytropes, we note finally that, for D > 2, Eq. (jl()2|) is undetermined for n = n 5 . Calculating 
the kinetic energy K = (D/2) J pd D r with Eqs. (|HU|) . (jHSJ), and using the Virial theorem 
(}9"T]) to obtain the potential energy, we find after simplification that 



(HI) 



A, 



T 



(4-£>) 



1 + 



or 



D(D 



D 



D+2 



a 



h + € 2 1 



2^2 



For D = 3, the integral is explicitly given by 

e 



(112) 
We note that for a 



9a{a 2 - 3) 3^ 



( a\ 
arctan —= . 

\V3j 



o (1 + f) 3 ^ 8(a> + 3)=> 
+oo, the energy diverges like A5 ~ (nV3/64)a for D = 3. 
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3.6 The minimum temperature and minimum energy 

The curve 77(a) presents an extremum at points a^ such that dr]/da{a^) = 0. Using Eqs. (|S7f). 
(|7T|) and (J72J), we find that this condition is equivalent to 

, . (D-2)n-D 

(113) % = ^ ^ =ii.. 

n — 1 

Therefore, the points where 77 is extremum are determined by the intersections between the 
solution curve in the (u,v) plane and the straight line defined by Eq. IjllHj) . The number 
of extrema depends on the value of D and n. It can be determined easily by a graphical 
construction using Figs. [IE] (an explicit construction is made in Fig. EDI see also [29J for D = 3). 
For D < 2, u s < for n > 1 and u s > D for n < 1 so that there is no extremum (case A). For 
2 < D < 10, there is no extremum for n <n^ (case A), there is one maximum for 77,3 < n < 
(case B) and there is an infinity of extrema for n > n§ (case C). They exhibit damped 
oscillations towards the value rj s corresponding to the singular solution (|47p. Asymptotically, 
the a>k follow a geometric progression a k ~ exp{2fc7r/ \J — A} (see for D = 3). For D > 10, 
there is no extremum for n < n 3 (case A), there is one maximum for n 3 < n < n 5 (case 5), 
there is an infinity of extrema for n 5 < n < n_ (case C) and there is no extremum for n > n_ 
(case D). This last case corresponds to an overdamped evolution towards the value r] s (in our 
mechanical analogy of Sec. 13. 3J) . For n = n^, using Eq. ([90)1 . we find that the extremum is 
located at «i = \J D(D + 2). For incomplete polytropes, the parameter rj is restricted by the 
inequalities (see Figs. EEl) 

(114) f] < uj n (case A), 

(115) rj < rj(ai) (cases B and C), 



(116) rj < Tj s (case D). 

These inequalities determine a maximum mass (for given T and R) or a minimum temperature 
(for given M and R) beyond which the system will either converge toward a complete polytrope 
with radius R* < R (if it is stable) or collapse. This dynamical evolution will be studied in 
Sec. E]for the nonlinear Smoluchowski-Poisson system. 

The curve A(a) presents an extremum at points a' k such that dA/da(a' k ) = 0. Using 
Eqs. f|102J) . (J7TJ) and ([72)1. we find that this condition is equivalent to 

2(n+ 1 - D)ul + {n + l)(n + 1 - D)u v + 2(D-n - 1)(D - l)u 

(117) +^D(D - 4)(n + l)u + ~D{D - 4)(n + l)u + ^D(D - 4) (2 - D)(n + 1) = 0. 

We can check that the point (it s , f s ) is a solution of this equation. On the other hand, vq = D — 2 
for Mo = and v ~ —2uo/(n + 1) for u — > ±oo. Finally, t> o — > oo for u — > where 

, s D(4-D) 
118 it* - 



2(n + l - D) 
More precisely, for u — > we have 

. . . . D(D-4)(D-2) . w . . , 

(119) (wo - u*)v ~ — — — — — ( n - n 3 / 2 )(n - ri 5 ) = A a (n). 
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Figure 5: Mass-central density profiles for polytropic configurations in a space of dimension 
2 < D < fO (specifically D = 3). A mass peak appears for the first time for the index n^. For 
n > ?25 the profile displays an infinity of peaks. For n < 1, 77 is a decreasing function of a. For 
D < 2, the 77(a) curves are monotonic. 




Figure 6: Mass-central density profiles for polytropic configurations in a space of dimension 
D > 10 (specifically D = 15). For n 5 < n < n_ (specifically n = 1.321) the profile displays 
an infinity of peaks. For n > n_ (specifically n = 2.5) the function 77(a) — > 77^ for a — > +00 
without oscillating. 
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Figure 7: Graphical construction determining the extrema of A (a) for D < 2 (specifically 
D = 1). The solid lines correspond to the solution curves and the dashed lines to the curves 
defined by Eq. (|117|) . The curves are labeled by the value of the polytropic index n. The 
vertical lines correspond to the asymptote u = u*. For D < 2, there is no intersection so that 
A(a) has no extremum. 
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Figure 9: Same as Fig. for 4 < D < 2(1 + \/2) (specifically D = 4.5). The geometrical 
construction changes for n = n 3 / 2 , n = and n = D — 1. Typical cases are represented. The 
indices label both the solid curve and the closest broken curve. 




Figure 10: Same as Fig. Elfor 2(1 + V2) < D < 10 (specifically D = 6). The geometrical 
construction changes for n = n 5 , n = TI3/2 and n = D — 1. Typical cases are represented. The 
indices label both the solid curve and the closest broken curve. 
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Figure 11: Evolution of the energy along the series of equilibria (parameterized by a) for 
2 < D < 4 (specifically D = 3). For n < n 5 , the curve has no extremum. For D < 2, the A(a) 
curves are monotonic. 




Figure 12: Evolution of the energy along the series of equilibria (parameterized by a) for 
4 < D < 10 (specifically D = 4.5). For n < n 5 , the curve has one maximum. For D > 10, the 
A (a) curves are similar to the t](a) curves in Fig. |H1 
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Figure 13: Generalized caloric curve for D < 2 (specifically D — 1). Note that according to 
Eq. ()98|). the potential energy is necessarily positive for D < 2, so the region A > is forbidden. 
We have plotted in dashed line the branch of complete polytropes with < R defined by Eq. 
flEED - 




Figure 14: Generalized caloric curve for D = 2. We have plotted in dashed line the branch of 
complete polytropes with < R defined by Eq. (JTTUJ). 
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Figure 15: Generalized caloric curve for 2 < D < 4 (specifically .D = 3). For 77,3 < n < 77,5, 
the inverse temperature presents a maximum but not the energy. For n > n 3 , there exists a 
region of negative (generalized) specific heats C = dE/dT < in the microcanonical ensemble. 
We have plotted in dashed line the branch of complete polytropes with R* < R defined by Eq. 
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Figure 16: Generalized caloric curve for 4 < D < 10 (specifically D = 4.5). For n < 
(specifically 7^3 = 1.8), the energy presents a minimum but not the inverse temperature. For 
ris < n < 723/2 (specifically 77,3/2 = 2.25), both the energy and the temperature present a 
minimum and the caloric curve 77(A) rotates anti-clockwise. This implies that equilibrium 
states with positive as well as negative specific heats C = dE/dT are stable in the canonical 
ensemble. This "thermo dynamical anomaly" arises because, as discussed in Sec. I3.1( stellar 
polytropes with n < n 3 / 2 are unphysical (the temperature is negative). For n = n 3 / 2 (white 
dwarfs), the curve makes an angular point (see Appendix lUj) . 
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Figure 17: Continuation of Fig. For n 3 / 2 < n < n 5 (specifically n 5 = 2.6) both the energy 
and the temperature present an extremum and the curve rotates clockwise (the curve makes a 
"loop"). The region of negative specific heats is now unstable in the canonical ensemble, as it 
should. For n > n 5 (specifically n 5 = 2.6), the energy and temperature present an infinity of 
extrema. 




Figure 18: Generalized caloric curve for D > 2(1 + y/2) (specifically D = 5.1) and n = n 3 / 2 - 
For this particular index, the curve presents an infinity (because n 3 / 2 > n 5 ) of angular points 
towards the singular sphere (see Appendix [CJ). 
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Figure 19: This figure summarizes the structure of the caloric curve as a function of the 
dimension D and the polytropic index n. The symbols are defined in the text. If we limit 
ourselves to n > D/2 (physical stellar polytropes), the region AB' showing a "thermodynamical 
anomaly" is not accessible. 



The two roots of A2(n) are n = 113/2 and n = n*,. They coincide at the particular dimension 
D — 2(1 + y/2). We note also that u* = for D = 4. For D > 2 and n = D - 1, Eq. (fTTTI) 
reduces to u + v = D — 2. 

The points where A(a) is extremum are determined by the intersections between the solution 
curve in the (u, v) plane and the curve defined by Eq. (|117j) . The number of extrema can thus 
be determined by a graphical construction using Figs. HHU This graphical construction depends 
on the values of D and n and the different cases are shown in Figs. I7IITU1 For D < 2, there is 
no extremum (case A'). For 2 < D < 4, there is no extremum for n < (case A') and there is 
an infinity of extrema for n > n 5 (case C). They exhibit damped oscillations toward the value 
A s corresponding to the singular solution (J47)) . For 4 < D < 10, there is one maximum for 
n < n§ (case B') and there is an infinity of extrema for n > n 5 (case C). For D > 10, there is 
one maximum for n < n§ (case B'), there is an infinity of extrema for n 5 < n < (case C), 
and there is no extremum for n > n_ (case D'). This last case corresponds to an overdamped 
evolution towards the value A s (in our mechanical analogy of Sec. 13. 3|) . The appearance of 
a maximum for n < n 5 when D > 4 was a surprise in view of preceding analysis for D = 3 
[2T| |2H1 12H] • The parameter A is restricted by the inequalities (see Figs. [TT] and ["HZ]) 



(120) 



A < A, 



(case A'), 



(121) 



A < A(ai) 



(cases B' and C) 



(122) 



A < A. e 



(case D'). 



These inequalities determine a minimum energy (for given M and R) below which the system 
will either converge toward a complete polytrope with radius -R* < R (if it is stable) or collapse. 
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In Figs. TT3WTSI we have plotted the generalized caloric curve A — rj, giving the inverse 
temperature as a function of the energy, for different dimensions D and polytropic index n. 
This extends the results of our previous analysis in D = 3 • The number of turning points 
as a function of D and n is recapitulated in Fig. 



3.7 The generalized thermodynamical stability 



We now come to the generalized thermodynamical stability problem. We shall say that a 
polytrope is stable if it corresponds to a maximum of Tsallis entropy (free energy) at fixed 
mass and energy (temperature) in the microcanonical (canonical) ensemble. The stability 
analysis has already been performed for D = 3 j23 12H I2H EH] and is extended here to a space 
of arbitrary dimension. This stability analysis can be used to settle either (i) the dynamical 
stability of stellar and gaseous polytropes (see Sec. |2J) or (ii) the generalized thermodynamical 
stability of self-gravitating Langevin particles (see Sec. HJ). 

We start by the canonical ensemble which is simpler in a first approach. A polytropic 
distribution is a local minimum of free energy at fixed mass and temperature if, and only if, 
the second order variations (see Appendix ID|) 



(123) 



5p5$d D r, 



are positive for any perturbation Sp that conserves mass, i.e. 



(124) 



8p d 



0. 



This is the condition of (generalized) thermodynamical stability in the canonical ensemble. 
Introducing the function q{r) by the relation 



(125) 



Sp 



1 



dq 



Spr 1 dr ' 



and integrating by parts, we can put the second order variations of free energy in the quadratic 
form 



(126) 



6 2 F 



R 



drq 



dr KSnr 1 dr 



G 



q- 



The second order variations of free energy can be negative (implying instability) only if the 
differential operator which occurs in the integral has positive eigenvalues. We need therefore 
to consider the eigenvalue problem 



(127) 



d ( p"< 



d 



dr \Srjr 



D-l 



dr 



+ 



G 



■ D-l 



q x (r) = Xqx(r). 



with q\ (0) = qx(R) = in order to satisfy the conservation of mass. If all the eigenvalues A are 
negative, the polytrope is a minimum of free energy. If at least one eigenvalue is positive, the 
polytrope is an unstable saddle point. The point of marginal stability in the series of equilibria 
is determined by the condition that the largest eigenvalue is equal to zero (A = 0). We thus 
have to solve the differential equation 



12? 



d_( _p 



dF 



dr KSdt 1 dr 



+ 



GF 

r D-l 



0. 
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with .F(O) = F(R) = 0. The same eigenvalue equation is obtained by studying the linear 
stability of the Euler- Jeans equation [23 EH] • Introducing the dimensionless variables defined 
previously, we can rewrite this equation in the form 

, N d f0 l - n dF\ nF 

(129) -A—— )+775TT = 0, 

with F(0) = F(a) = 0. If 



< 130 > ^TAi^TdU^ 

denotes the differential operator occurring in Eq. ()129j) . we can check by using the Emden 
Eq. (gij) that 

(131) C{£ D - 1 &) = {n-l)&, C(£ D 6 n ) = [(2 — D)n + D]6'. 

Therefore, the general solution of Eq. ()129j) satisfying the boundary conditions at £ = is 

(132) F(() = c 



er + (g - 2 »" - y- v 

n — 1 



Using Eq. ()132|) and introducing the Milne variables (|70j) . the condition F(a) = can be 
written 

(g — 2)n — D 

(133) m = = u s . 

n — 1 

This relation determines the points at which a new eigenvalue becomes positive (A = + ). 
Comparing with Eq. ()113|) . we see that a mode of stability is lost each time that rj is extremum 
in the series of equilibria, in agreement with the turning point criterion of Katz jJT] in the 
canonical ensemble. When the curve rj(a) is monotonic (cases A and D), the system is always 
stable because it is stable at low density contrasts (a —>■ 0) and no change of stability occurs 
afterward. When the curve r](a) presents extrema (cases B and C), the series of equilibria 
becomes unstable at the point of minimum temperature (or maximum mass) ot\. In Fig. E3 
this corresponds to a point of infinite specific heat C = dE/dT — > oo, just before entering the 
region of negative specific heats C < 0. When the curve 77(a) presents several extrema (case 
C), secondary modes of instability appear at values a 2 , a 3 ,... (see in D — 3). We note that 
complete polytropes (with n < if D > 2) are stable in the canonical ensemble if D < 2 and 
if (D > 2, n < 723). They are unstable otherwise. In the thermodynamical analogy developed in 
|29| I26|. this is a condition of nonlinear dynamical stability for gaseous polytropes with respect 
to the Jeans-Euler equations. In particular, the self-gravitating Fermi gas at zero temperature 
(a classical white dwarf star) is dynamically stable if n 3 / 2 < n 3, i-e. D < 4, and unstable 
otherwise. 

According to Eq. (|125|) . the perturbation profile that triggers a mode of instability at the 
critical point A = is given by 



(134) 



5p 1 dF 



Po S D £ D -i d£ ' 

where F(£) is given by Eq. ()132j) . Introducing the Milne variables (|70j). we get 
(135) — = — -{y s -v). 
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Figure 20: Location of the turning points of temperature in the (u, v) plane for systems with 
dimension 2 < D < 10 (specifically D = 3). The line u = u s determines the extrema of r\ and 
the line v — v s determines the nodes of the density profiles that trigger the instabilities in the 
canonical ensemble. 



The density perturbation 5p becomes zero at point(s) £j such that = v s . The number of 
nodes is therefore given by the number of intersections between the solution curve in the (u, v) 
plane and the line v = v s . It can be determined by straightforward graphical constructions in 
the Milne plane, using Figs. HHU (see, e.g., Fig. EH]) . When the solution curve is monotonic 
(case B), the density perturbation profile has only one node. In particular, for n = n^, the 
perturbation profile at the point of marginal stability is given by 



(136) 



5p _ (D + 2)ci D(D- 2) - £ 2 
7 ~ 2S D D{D-2)+e' 



It vanishes for £(i) = a/ D(D — 2). When the solution curve forms a spiral (case C), the density 
perturbation 5p corresponding to the k-th mode of instability has k zeros £i, £ 2) •■■> ^fc < a k- In 
particular, the first mode of instability has only one node. For high modes of instability, the 
zeros asymptotically follow a geometric progression with ratio 1 : exp{27r/v^-A} (see Ref. |2E] 
for D = 3). 

In the microcanonical ensemble, a polytrope is a maximum of entropy at fixed mass and 
energy if, and only if, the second order variations (see Appendix [01 



(137) 



2n 1 
D(2n-D) J pd D r 



p 2 2 



$ + ^)5pd D r 
2 pj 



are negative for any variation Sp that conserves mass to first order (the conservation of energy 
has already been taken into account in obtaining Eq. (|137|0 . Now, using Eq. (|125j) and 
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integrating by parts, the second variations of entropy can be put in a quadratic form 

-R r R 



(138) 
with 



5 2 S 



o Jo 



drdr'q(r)K(r, r')q(r') 



K(r } r') 



2n 



D(2n -D) J pd D r 



(r) ( $ + jK^p^- 1 



(r') 



(139) 



+~5(r-r') 



" , d ( p 7 " 2 d 
A7— 



dr \Srjr D 1 



G 



,D-1 



The problem of stability can therefore be reduced to the study of the eigenvalue equation 

-R 



(140) 



dr'K(r,r')q x (r') = Xq x {r), 



with q\(0) = q\(R) = 0. The point of marginal stability (A = 0) will be determined by solving 
the differential equation 



(141) 
with 
(142) 



T d( p 7 " 2 dF 
K 7— 



dr \ Sdt 1 dr 



+ 



GF 

~.D-1 



K 1 2 /o DM/ 7-2 d P 



V 



In arriving at this expression, we have used the relation 



(143) 



In dr 



which results from the condition of hydrostatic equilibrium (jl3j) with the polytropic equation 
of state (J3n|l . Introducing the dimensionless variables defined previously, Eqs. (j!41j) and ()142j) 
can be rewritten 



(144) 

with 

(145) 



fl 1 " 71 dF 



d ( 

~dt\t D - x dt) 1 e- 1 



y = —(n + l)(2n- D) ra Jo - „ 1 — . 



and -F(O) = -F(ci) = 0. Using the identities (|131|) . we can check that the general solution of 
Eq. ()144|) satisfying the boundary conditions for £ = and £ = a is 



(146) 



no 



D-ln/\ 



(n - l)uo + D — (D — 2)n 



The point of marginal stability is then obtained by substituting the solution ()14(ij) in Eq. (|145j) . 
Using the identities (see Appendix |EJ) 



(147) 



D-1/qi\2 



di 



a D 6'(a) 2 



D + 2 - (D - 2)n 
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(149) / fl W = „ ° C "' ( ° )2 ; f-J - (J - 2)3! + ~ 2 > 



D + 2 - (£> - 2)n V v v 

which result from simple integrations by parts and from the properties of the Lane-Emden 
equation (JUJ, it is found that the point of marginal stability is determined by the condition 
(|117|) . Therefore, the series of equilibria becomes unstable at the point of minimum energy 
in agreement with the turning point criterion of Katz [H] in the microcanonical ensemble. 
When the curve A(a) is monotonic (cases A' and D f ), the system is always stable. When the 
curve A(a) presents extrema (cases B' and C), the series of equilibria becomes unstable at the 
point of minimum energy a[. In Fig. E3 this corresponds to the point where the specific heat 
C = dE/dT = 0, passing from negative to positive values. Note that the branch of negative 
specific heats between CE and MCE is stable in the microcanonical ensemble although it is 
unstable in the canonical ensemble. When the curve A(ct) presents several extrema (case C), 
secondary modes of instability appear at values a' 2 , ot' 3 ,... We note that complete polytropes 
(with n < ris if D > 2) are stable in the microcanonical ensemble if D < 4 and unstable if 
D > 4. Owing to the thermodynamical analogy, this is a condition of nonlinear dynamical 
stability for stellar polytropes with respect to the Vlasov equation [SHj. The difference between 
the dynamical stability of gaseous polytropes (n < n 3 ) and stellar polytropes (n < n 5 ) was 
related in |2T?| EE] to a situation of ensemble inequivalence (and the existence of a negative 
specific heat region) in thermodynamics. Since the caloric curve is monotonic in D = 2, we also 
conclude that polytropic vortices [20] are always nonlinearly dynamically stable with respect 
to the 2D Euler equation. 

According to Eqs. (J 134)) and (|146j) . the perturbation profile that triggers a mode of insta- 
bility at the critical point A = is given by 

(!50) — = — — jjr — -r-(D-nv-Ua), 

p Sd [n — 1)uq + D — [D — 2)n 

where we have used the Emden equation (}45|) and introduced the Milne variables (|70|) . The 
number of nodes in the perturbation profile can be determined by a graphical construction 
similar to the one described in Refs. [12113 f° r n = 00 (isothermal case). 



4 Self-gravitating Langevin particles 
4.1 The nonlinear Smoluchowski-Poisson system 

Let us consider a system of self-gravitating Brownian particles described by the stochastic 
equations (i = 1, N) 

(151) ^ = Vi ' ^ = -V$ l -^v,, + ^/2T|R l (t), 

where $j = $(r^,t) is the self-consistent gravitational potential created by the particles, — £vj 
is a friction force originating from the presence of an inert medium and Rj(t) is a white noise 
satisfying (Rj(i)) = and (Ri ta (t)Rj,b(t')) = 5ijS a bS(t—t'). To simplify the problem, we consider 
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the high friction limit £ — > +00, where £ is the friction coefficient [2]. This regime is achieved 
for times t ^> In the mean-field approximation, the evolution of the density of particles is 
governed by the Smoluchowski equation [TU] 



V(Dp) + ipV$ 



coupled to the Newton- Poisson equation (jlj). In the usual case [21 HUE], the diffusion coefficient 
D is constant and the condition that the Boltzmann distribution p ~ e - *^ is a stationary 
solution of Eq. (J152j) is ensured by the Einstein relation £D = T. It can be shown |2j that the 
SP system decreases the free energy F = E — TS constructed with the Boltzmann entropy. 
Hence, the equilibrium state minimizes 

(153) F = ~ J p$d D r + T J p In pd D r , 

at fixed M and T. 

We want to consider here a more general situation in which the diffusion coefficient D 
depends on the density p while the drift coefficient £ is still constant. In the absence of drift, 
this would lead to a situation of anomalous diffusion. In the presence of drift, a notion of 
generalized thermodynamics emerges [20J. Indeed, writing the diffusion coefficient in the form 
D(p) = jp(p)/p, we obtain the generalized Smoluchowski equation 



i(Vp + pV$) 

In Ref. [20J, it is shown that generalized Smoluchowski equations of this type satisfy a form of 
canonical if-theorem. The Lyapunov functional decreasing monotonically with time is 



(155) F = J p J ^±dp'd u r + -J p§d u r. 

This can be interpreted as a free energy associated with a generalized entropy functional (see 
[2"U] for more details). The equilibrium state minimizes F at fixed M. In the present context, 
this is a condition of generalized thermodynamical stability in canonical ensemble. 

The generalized Smoluchowski equation (|154j) can be obtained by combining the ordinary 
Fokker-Planck equation with a Langevin equation of the form [2*U] : 

where R(t) is a white noise. When p(p) is a power-law, Eq. ()156|) reduces to the stochas- 
tic equations studied by Borland [21] in connexion with Tsallis thermodynamics. Since the 
function in front of R(t) depends on r, the last term in Eq. (|156J) can be interpreted as a 
multiplicative noise. Note that the noise depends on r through the density p(r). Kaniadakis 
j22] has also introduced generalized Fokker-Planck equation arising from a modified form of 
transition probabilities. In these works, the Langevin particles evolve in an external potential. 
The case of Langevin particles in interaction was considered by Chavanis [20] • He introduced 
a generalized Fokker-Planck equation (see in particular Eq. (81) of 20 ) valid for an arbitrary 
equation of state p = p(p), or diffusion coefficient D(p), and for an arbitrary binary potential 
of interaction u(r — r'). This equation has been studied recently in PJ |2l El Ej for an isothermal 



34 



equation of state p = pT (constant D) and a gravitational interaction. This study has been 
extended in [13] to a Fermi-Dirac equation of state. In this paper, we consider the case where 
the function D(p) is a power law and write £D(p) = Kp 1 ^ 1 . This corresponds to a polytropic 
equation of state p = Kp 1 . Then, Eq. (j 152)1 can be rewritten 



(1B7) | = V 



^(iTVp 7 + pV$) 



For the nonlinear Smoluchowski-Poisson system, the Lyapunov functional decreasing monoton- 
ically with time is j2Dl 

(158) F = — [ (p 7 - p)d D r + - I p$d D r. 

7 - 1 J 2 J 

It can be interpreted as a free energy associated with Tsallis entropy. In this context, the 
polytropic index 7 plays the role of the g-parameter and the polytropic constant K plays the 
role of the temperature (see Ref. [20] for details and subtleties). Therefore, keeping K fixed 
corresponds to a canonical situation [29 . For 7 — >• 1, we recover the Boltzmann free energy 
studied in [21 EJ. For 7 = 5/3, i.e. n = 3/2, Eq. ()157|) describes self-gravitating Brownian 
fermions at T = (in D = 3) g3|. 

The nonlinear Smoluchowski equation can be obtained from a variational principle, called 
Maximum Entropy Production Principle, by maximizing the rate of free energy F for a fixed 
total mass M [HI2n]- It is straightforward to check that the rate of free energy dissipation can 
be put in the form 

(159) F = — / ^-{KVp 1 + pV$) 2 d D r < 0. 

J pi 

For a stationary solution, F — 0, and we obtain a polytropic distribution which is a critical 
point of F at fixed M. Considering a small perturbation around equilibrium, we can establish 
the identity [Hj: 

(160) 5 2 F = 2X5 2 F < 0, 

where A is the growth rate of the perturbation defined such that 5p ~ e A *. This relation shows 
that a stationary solution of the nonlinear Smoluchowski-Poisson (NSP) system is dynamically 
stable for small perturbations (A < 0) if and only if it is a local minimum of free energy 
(5 2 F > 0). In addition, it is shown in Refs. [21120] that the eigenvalue problem determining the 
growth rate A of the perturbation is similar to the eigenvalue problem ()128j) associated with 
the second order variations of free energy (they coincide for marginal stability). This shows the 
equivalence between dynamical and generalized thermodynamical stability for self-gravitating 
Langevin particles exhibiting anomalous diffusion (this result was obtained independently by 
Shiino 24] in the specific context of Tsallis thermodynamics). In fact, our formalism is valid for 
more general functionals than the Boltzmann or the Tsallis entropies [201- These functionals 
()155|1 arise when the diffusion coefficient is of the general form D(p), not necessarily a power 
law. We note that the NSP system satisfies a form of Virial theorem [2U]. For D ^ 2, it reads 

(161) ^ = 2K + (D — 2)W - Dp b V, 
where 

(162) I=\ pr 2 d D r, 
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is the moment of inertia and pb the pressure on the box (assumed uniform). For D = 2, the 
term (D — 2)W is replaced by —GM 2 /2. For a stationary solution dl/dt = 0, we recover the 
Virial theorem (|9ip. 

Our model of self-gravitating Brownian (or Langevin) particles has no clear astrophysical 
applications, so it must be regarded essentially as a toy model of gravitational dynamics. It may 
find application for the formation of planetesimals in the solar nebula since the dust particles 
experience a friction with the gas and a noise due to small-scale turbulence |44j . However, 
even in this context, the model has to be refined so as to take into account the attraction of 
the Sun and the rotation of the disk. Anyway, the self-gravitating Brownian (or Langevin) 
gas model is well-posed mathematically, and it possesses a rigorous thermodynamical structure 
corresponding to the canonical ensemble. Therefore, it can be used as a simple model to 
illustrate some aspects of the thermodynamics of self-gravitating systems. Since it minimizes 
F = W[p] at fixed mass, it can also be used as a powerful numerical algorithm to construct 
nonlinearly dynamically stable stationary solutions of the Euler- Jeans equations (see Sec. I2.2JI . 
Coincidentally, the SP system also provides a simple model for the chemotactic aggregation 
of bacterial populations |TT] . The name chemotaxis refers to the motion of organisms induced 
by chemical signals. In some cases, the biological organisms secrete a substance that has an 
attractive effect on the organisms themselves. This is the case for the bacteria Escherichia coli. 
In the simplest model, the bacteria have a diffusive motion and they also move systematically 
along the gradient of concentration of the chemical they secrete. Since the density $(r, t) 
of the secreted substance is induced by the particles themselves, the drift is directed toward 
the region of higher density. This attraction triggers a self-accelerating process until a point 
at which aggregation takes place. If we assume in a first step that $(r, t) is related to the 
bacterial density p(r, t) by a Poisson equation, this phenomenon can be modeled by the SP 
system. Now, it has been observed in many occasions in biology that the diffusion of the 
particles is anomalous [llj. This is a physical motivation to study the NSP system in which 
the diffusion coefficient is a power law of the density. In the following section, we show that the 
NSP system admits self-similar solutions describing the collapse of the self-gravitating Langevin 
gas or of the bacterial population. These theoretical results are confirmed in Sec. El where we 
numerically solve the NSP system. 



4.2 Self-similar solutions of the nonlinear Smoluchowski-Poisson sys- 
tem 

From now on, we set M = R = G = £ = 1 without loss of generality. The equations of the 
problem become 

(163) % = V(iTVp 7 + pV$), 



(164) A$ = S DP , 
with boundary conditions 

(165) ^(0,t)=0, *(1) = ^, ^(1)+ P (1) = 0, 

for D^2. For D = 2, we take $(1) = on the boundary. We restrict ourselves to spherically 
symmetric solutions. Integrating Eq. (|164j) once, we can rewrite the NSP system in the form 
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of a single integrodifferential equation 

1 d 



dp 
dt 



(166) 

where we have set 
(167) 
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The quantity 9 can be seen as a sort of generalized temperature (sometimes called a polytropic 
temperature jSH]) and it reduces to the ordinary temperature T for n — > +oo. We note that the 
proper description of a gas of Langevin particles in interactions is the canonical ensemble where 
is fixed. However, we can formally set up a microcanonical description of self-gravitating 
Langevin particles by letting the temperature 0(i) depend on time so as to conserve the total 
energy: 



(168) 



^=f=^/e(«)^r + i/^. 



The two situations have been considered in the case of self-gravitating Brownian particles 
(n^+oo) in [UI21I3I1I. 

The NSP system is also equivalent to a single differential equation 
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~~dT 



1 dM 



,D-1 
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d 2 M D-l dM 
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r dr 
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M dM 



,D-1 



dr 



M{r,t) = / p{r f ,t)S D r' D - 1 dr', 



(169) 

for the quantity 
(170) 

which represents the mass contained within the sphere of radius r. The appropriate boundary 
conditions are 

(171) M(0,t) = 0, M(l,t) = l. 

It is also convenient to introduce the function s(r,t) = M(r,t)/r D satisfying 



(172) 
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di 



B [ r |£ + jDs 
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r dr 



+ f + Ds ).s. 

dr 



For n — > +oo, these equations reduce to those studied in Refs. pOElll] in the isothermal case. 
We look for self-similar solutions of the form 



(173) 



p(r,t) = p (t)f 



ro(t) 
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1-1 In 



1/2 



The radius r defined by the foregoing equation provides a typical value of the core radius of 
an incomplete polytrope (with n > n^). It reduces to the King's radius ^H] as n -> +oo. In 
terms of the mass profile, we have 



(174) 



M{r,t) = M {t)g 



ro(t)>/' 



with 



M Q (t) = p Q r L Q 



D 
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and 

(175) g(x) = / f(x')S D x' D ~ l dx'. 

Jo 

In terms of the function s, we have 

(176) s(r,t) = p (t)s(^—\ with S(x" 



r o(t) J ' x D 

Substituting the ansatz (|176|) into Eq. ()172|) . and using the definition of r in Eq. ()173|) . we 
find that 

(177) ^ S - ^^xS' = plixS' + DS) l ' n ( 3" + ^±S') + p 2 JxS' + D3)S, 
at r at ' \ x J 

where we have set x = v/vq. We now assume that there exists a such that 

(178) p ~ r a . 
Inserting this relation into Eq. ()177j) . we find 



(179) 



^r(s+ ^xS'^j = pi (xS' + DS) 1/n (s" + + (xS' + DS)S 



which implies that (1 / p^) (dp / dt) is a constant that we arbitrarily set equal to a. This leads 
to 

(180) Po (t) = -{tcou-ty 1 , 

a 

so that the central density becomes infinite in a finite time t co u. The scaling equation now reads 

(181) aS + xS' = (xS' + DS) 1,n (V + ^Is^j + (xS 1 + DS)S. 
For x — > +oo, we have asymptotically 

(182) S(x) ~ x~ a , g{x) ~ x D ~ a , f(x) ~ 

In the canonical ensemble where 6 is constant, Eq. ()173j) and Eq. ()178|) lead to a = a n , 
with 

(183) a 



n -. • 

71—1 

Note that for n — > oo, we recover the result of = 2. Equation ()182j) implies that for 

large x, p ~ (Z) — a)5 > 0, which enforces a < D (this also guarantees that the mass of 
the power-law profile p = Cr~ a at t — t co u is finite). The limit value a n = D corresponds 
to n = n 3 . Therefore, there is no scaling solution for n < n 3 . This is consistent with our 
finding that the collapse occurs only for n > n 3 . For n < n 3 and rj > u n , the system converges 
toward a complete polytrope with radius R* < 1 which is stable (see Sec. 13.71 and Fig. [TBI . 
For n 3 < n < and rj > r/(ai), we can formally construct a complete polytrope with radius 
i?* < 1, but this structure is unstable (Sec. 13. 7j) so that the system undergoes gravitational 
collapse. 
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In the microcanonical ensemble, the value of a > a n cannot be obtained by dimensional 
analysis. It will be selected by the dynamics. In the case n — > +00, we have found in [3| 
that Eq. (J181J) for the scaling profile has physical solutions only for 2 < a < a max (-D) (with 
tt m ax(-D = 3) = 2.209733...). For arbitrary n, such a a max (-D, n) > a n also exists (see Sec. EJ). It 
is easy to see that the maximum value for a leads to the maximum divergence of temperature 
and entropy. Therefore, it is natural to expect that the value a max will be selected by the 
dynamics except if some kinetic constraints forbid this natural evolution (see below). In fact, 
as already noted in [2j, a value of a > a n poses problem with respect to the conservation of 
energy. We recall (and generalize) the argument below. According to Eqs. (j!73|) . (j!78|) and 
()180j) . during collapse, the temperature behaves as 

(184) ~ pl- 1/n - 2/a ~ (t coll - *)-<>/«»-*/«>, 
and the kinetic energy ()34j) behaves as 

(185) K~Q f p\r,t)r D - 1 dr ~ 9(p ro) 7 x I r ' 1 '^ dr ~ 6 / r ' 1 ^ dr. 

J J VQ J TQ 

First consider the case n > n 5 . If a < D/7 (which is the case in practice since a n < -D/7 
implies n > n 5 ), the integral is finite and the kinetic energy behaves as K ~ 0. Therefore, it 
diverges at t co u for any a > a n . On the other hand, the scaling contribution to the potential 
energy behaves as 

(186) W ~ f dr ~ (p r£) 2 x f r D+1 ~ 2a dr ~ f r D+1 ~ 2a dr. 

If a < (D + 2)/2 (which is the case in practice since a n < (D + 2)/2 implies n > 715), the scaling 
part of W remains finite at t co u. Energy conservation would then imply that a = a n . In a 
first series of numerical experiments reaching moderately high values of the central density [2] , 
we measured (by different methods) a scaling exponent a ~ 2.2 > = 2 (for the isothermal 
case n = 00). Combined with the fact that the Smoluchowski-Poisson system must lead to a 
diverging entropy, we argued that a max is selected by the dynamics (while being careful not 
to rigorously reject the possibility that a = = 2). Then, in order to account for energy 
conservation, we proposed a heuristic scenario showing how subscaling contributions could 
lead to the divergence of potential energy. In fact, the numerical simulations were not really 
conclusive in showing the divergence of temperature, as the expected exponent is very small 
(2/a n — 2/a max = 0.09491..., for D = 3 and n — ► +00). Recently, we have conducted a new 
series of numerical simulations allowing to achieve much higher values of density (see Sec. |3J). 
These simulations tend to favor a value of a = a n leading to a finite value of the temperature 
at t co u. However, the convergence to the value a n is obtained for values of t very close to t co u 
and, at intermediate times for which the temperature has not reached its asymptotic value, the 
numerical scaling function tends to display an effective exponent between a n and a max - This 
situation is reminiscent of the (D — 2,n — +00) case studied in 0, although the situation 
is not exactly equivalent. Numerical simulations of the (D = 3,n = +00) case have been 
conducted independently by P and favor also a value of a = a n . Note that there is no rigorous 
result proving that a = a n in the microcanonical situation, so this point remains an open 
mathematical problem. 

For n < 715, the kinetic and potential energies diverge in a consistent way as 

(187) K ~ 6 I r - 1 -^ dr ^ Qr^ a ~ pl~ (D+2)/a , 

Jr 

(188) W ~ !\- D+l - 2a dr~r° +2 - 2a ~ p l-^l a . 

J TQ 
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However, in the microcanonical ensemble, the system is expected to reach a self-confined poly- 
trope for n < n$ and A > A n since it is stable (see Sec. 13.71 and Fig. [T5]) . Probably, the choice 
of evolution will depend on a notion of basin of attraction as in [2|. 

We now focus on self-gravitating Brownian particles (n = oo). In case of collapse, the 
previous discussion shows that the system has the desire to achieve a value of a > 2 leading to a 
divergence of temperature and entropy. This is indeed a natural evolution in a thermodynamical 
sense. This is also consistent with the notion of gravothermal catastrophe introduced in the 
context of globular clusters [121 EE 113 EE] ■ However, the energy constraint (jl68J) seems to 
prevent this natural evolution (the divergence of entropy occurs in the post-collapse regime 0]). 
This is related to the assumption that the temperature is uniform although this assumption 
clearly breaks down during the late stage of the collapse. Therefore, we expect that if the 
temperature is not constrained to remain uniform, the system will select a value of a > 2 as 
in other models of microcanonical gravitational collapse ^OJEZIEE!- Below, we give a heuristic 
hint so as how this can happen. We consider the Smoluchowski equation 

(189) |Uv[V(Tp)+pW], 

where T = T(r,t) is now position dependant. In any model where the temperature T(r,t) 
satisfies a local conservation of energy, we expect the following scaling 



(190) T(r,t) =T (t)6[ —-), with 9(x) ~ x 2 ~ a , when x -» +oo. 

Such a scaling is indeed observed in the globular cluster model of [IB]. The decay exponent is 
obtained by using the definition of T ~ p rl and the fact that the temperature at distances 
r ^> r should be of order unity. The precise density and temperature profiles and the value of 
a depend on the model considered for the energy transport equation. It is not our purpose to 
discuss a precise model in the present paper and we postpone this study for a future work [30 . 
However, we can give an analytical argument showing why a > 2 should now be selected in a 
unique way. Equation (|190|) shows that temperature scales as the potential <3> since 



(191) $(r,t) =$(1) - / s{r',t)r'dr' ^ <t>(l) - p r 2 S{x)xdx^-T S{x)xdx. 

J t Jr/ro Jr/ro 

The scaling function for $(r, t) is then 

p+oo 

(192) <f>(x) = - S(x')x'dx'~x 2 - a , when x -> +oo. 



Equations (|190p . ()191[) and ()192|) imply that both the kinetic and potential energies remain 
bounded for all times t < t co u, at least for a < (D + 2)/2, even if the central temperature 
T (t) diverges. Indeed, the temperature increases in the core but the core mass goes to zero so 
that the kinetic energy of the core ~ M T goes to zero. On the other hand, the temperature 
remains of order unity in the halo leading to a finite kinetic energy in the halo. This "core- 
halo" structure for the temperature is more satisfactory than a model in which the temperature 
is uniform everywhere, even in the collapse phase. Before introducing a precise equation for 
T(r,t) [30J, we make the reasonable claim that the temperature and the potential energy are 
simply proportional in the core region as they exhibit a similar scaling relation. Defining T (t) 
such that #(0) = 1, we end up with the hypothesis 

(193) 6{X) = W) = ^C ° S{x ' )x ' dx '> 
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with 



D _ D 

(194) ~ ~W) ~ f +c °S(x')x>dx>' 
Using Eq. . we find that the scaling equation for S(x) is now 

(195) aS + xS' = ^4>(x) (s" + ^j-^'J + M - S( x> s' + DS). 

For a given A, this equation is an eigenvalue problem in a. In the limit of large dimension and 
proceeding exactly along the line of [3], it can be seen that S ~ 0(D _1 ) and A ~ 0(D°). Using 
the method of S(x) and a can be computed easily up to order 0(D -1 ), as a function of A 
and z = DS(0)/2. Now imposing the constraint of Eq. (j!94j) . this selects a unique value for a. 
After straightforward calculations, we find a simple parametrization of A and a as a function 
c£z = DS(0)/2 > 2 : 

(196) a -2 = - (1 -2Z- 1 ) + 0(ZT 2 ), 

(197) A = 2 (1 -z- 1 ) (1 -2Z- 1 ) + 0(ZT 1 ), 

which can be recast in the form 

(198) a -2 = A (VrT4A- l) +0(D- 2 ), 

(199) 5(0) = + 0(D- 2 ). 

Equation (|197|) implies that for A > 2 (up to order 0(D -1 )), there is no solution to the scaling 
equation, and that for A < 2, there is a unique a corresponding to a physical solution. In 
general, the actual value of A will be selected dynamically by an additional evolution equation 
for the temperature profile. However, assuming A ~ 1 is natural (although this point needs to 
be confirmed), since it corresponds naively to a local energy conservation condition 

(200) E T (r )t ) ~ -i$(r,t). 
In that case, we obtain 

(201) a -2 = A (VE- l) +0(ZT 2 ). 

The above argument is a strong indication that if the uniform temperature constraint is aban- 
doned, a non trivial value for a > 2 will be selected. We conjecture that this eigenvalue will 
be close to a max ~ 2.21 as found in other models of microcanonical gravitational collapse with 
a non-uniform temperature [JUl El OH] • 



5 Numerical simulations 

In this section, we illustrate numerically some of the theoretical results presented in the previous 
sections, but we restrain ourselves to D = 3. 

We first consider the dynamics in the canonical ensemble (fixed 0). In Fig. we show 
the different steps of the formation of a self-confined polytrope of index n = 3/2<ri3 = 3 
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Figure 21: Evolution of the density profile p(r,t) for n = 3/2 and = 0.1 (corresponding to 
rj > u n ). The profile converges to a complete polytrope strictly confined inside the box (thick 
line). The dashed line is the initial uniform density profile. 

similar to a classical "white dwarf" star. In this range of n, the system always converges to 
an equilibrium state (see Fig. H5|) . If rj < u n the equilibrium state is confined by the box 
(incomplete polytrope) while for 77 > uj n the density vanishes at < R (complete polytrope). 
In Fig. |22 w e illustrate the collapse dynamics at low temperatures for n = 4 6 [n 3 , n 5 ] and 
n = oo > n 5 = 5. This is compared to the predicted scaling profiles. The convergence to 
scaling is slower for n = 4 (a = 8/3) than for n = oo (a = 2). This is expected since, in the 
former case, tq ~ p decreases more slowly as a is bigger. Thus, the scaling regime tq ^ 1 
is reached in a slower way. For instance, for comparable final densities of order 10 6 , and for the 
considered temperatures, we find that the minimum r obtained for n = 4 is roughly 4 times 
bigger than in the n = oo simulations. For n = oo and in the large D limit, we have shown in 
[11] that the scaling function S(x) takes the form 



where x is such that S(x ) = a/D. The quantities z = DS(0)/2 and a(z) have been exactly 
calculated in this limit. In the present case, and for a given n (yielding a n — 2n/(n — 1)), 
we compute 5(0), by assuming the above functional form. The parameter Xq, X\ and S(0) 
are numerically calculated by imposing the exact value for S"(0) extracted from Eq. ()181|) . as 
well as the two conditions that xq must satisfy (see |3j for more details). The comparison of 
this approximate theory with actual numerical data is satisfactory (see Fig. I23J). Note that we 
have been unable to develop a large n perturbation theory in the same spirit as the large D 
expansion scheme derived in j^j. 

As explained in the previous section, the situation in the microcanonical ensemble (where 
= Q(t) evolves with time in order to conserve energy) is less clear. Like in the case n = oo 
studied in the scaling equation admits a physical solution for any a n < a < a ma _ x (D,n). 
In Fig. |2U we plot a max (-D = 3,n), as well as the corresponding value of 5(0), as a function 
of n. As explained previously, it is doubtful that a scaling actually develops with a > a n 



(202) 
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Figure 22: For n = 4 (a 4 = 8/3) and = 0.1 (canonical description), we plot S(r, t)/S(0, t) 
as a function of r/r (t), where r (t) is defined by Eq. ()173|) . for different times corresponding 
to central densities in the range 2.10 2 ~ 4.10 5 (bottom data collapse). This is compared to 
the scaling function obtained by solving Eq. (|181|) numerically (dotted line). The same is 
plotted in the case n = oo (a^ = 2), for which the scaling profile is known analytically [3] : 
S(x)/S(0) = (1 + x 2 )~ l (upper data collapse). The two curves have been shifted for clarity. In 
the n = oo case, the asymptotic scaling profile (dotted line) is almost indistinguishable from 
the data collapse. Dashed lines have respective slopes —8/3 and —2. 
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Figure 23: We plot S(0) as a function of n (full line), and compare it to a simple theory 
explained in the text, which is inspired by the large D perturbation introduced in [Hj (dashed 
line). Note that 5(0) ~ 4 + C/n + 0(n~ 2 ) for large n with C ~ 19. 
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Figure 24: We plot a max (D = 3, n) as a function of n (top plot), as well as the associated value 
of I In S(0) (bottom plot). The horizontal dashed line represents the asymptotic value of a max 
for n — > +00, and we find a max (D = 3, n) ~ a max (-D — 3,n — 00) + C^/n + 0{n~ 2 ) for large 
n, with C3 ~ 2.7. For n < n 5 = 5, a max = D~ = 3~ (strictly speaking, the scaling solution 
associated to a = D = 3 does not exist below n 5 = 5). We observe that In S(0) ~ (n — 5) -1 / 4 
provides an excellent fit of S(0) for n e [5, 10]. We have also plotted a n = 2n/(n — 1) (thick 
line). The scaling equation (|18ip admits solutions for a n < a < a max . The two curves intersect 
when a n = D which corresponds to n — n.3. There is no scaling collapse solution for n < n^. 
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Figure 25: For n = +00 and E = —0.45 (microcanonical description), we plot S(r,t)/S(0,t) as 
a function of r/ro(t) where ro(i) ~ po(£) 1//o , for times corresponding to central densities in the 
range 2.10 2 ~ 4.10 5 (for comparison, our previous simulations [2] did not exceed po ~ 1000). 
We try both values = 2 (bottom plot) and a = a max = 2.209733... (top plot), and compare 
both data collapses to the associated scaling function (dotted lines) . The two curves have been 
shifted for clarity. The scaling associated to a max is clearly more convincing than that for a = 2, 
especially at large distances. However, our simulations also suggest that 0(t) ~ p (t) 1_2//<:imax 
does not diverge at t co u (see the insert where a line of slope 1 — 2/a max ps 0.09491... has been 
drawn as a guide to the eye), so that the asymptotic scaling should correspond to a = 2. This 
apparent "paradox" clearly shows that the convergence to the limit value a = a n is extremely 
slow, suggesting an intermediate pseudo-scaling regime with a n < a < a max . 



when the temperature is uniform. However, a pseudo-scaling should be observed with a ~ 
a max . In Fig. |2H we present new simulations (D = 3, n = 00) confirming that the observed 
scaling dynamics is better described by a = a max than by a = 2, in the time/density range 
achieved. Such a value of a implies that the temperature would diverge with a small exponent 
(Q(t) ~ p (t) 1 ~ 2/amax , with 1 - 2/a max = 0.09491...). However, in the range of accessible 
densities po = 2.10 -1 ~ 10 6 , numerical data tend to suggest that the temperature converges 
to a finite value with an infinite derivative ^(t co u) = +00 as t ^ t co u. This convergence of 
the temperature has been observed independently by |9J. Thus, we conclude that the system 
first develops an apparent scaling with a < a max , before slowly approaching the asymptotic 
scaling regime with a = 2. In any case, it is clear that if the a = 2 scaling is the relevant one, 
the scaling regime is approached much more slowly than in the canonical ensemble (compare 
Fig. EH1 and Fig. 122(1 . This is a new aspect of the inequivalence of statistical ensembles for 
self-gravitating systems. 
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6 Conclusion 



In this paper, we have discussed the structure and the stability of self-gravitating polytropic 
spheres by using a formalism of generalized thermodynamics [20J. This formalism allows us to 
present and organize the results in an original manner. What we mean by generalized ther- 
modynamics is the extension of the usual variational principle of ordinary thermodynamics 
(maximization of the Boltzmann entropy Sb at fixed mass M and energy E) to a larger class 
of functionals (playing the role of "generalized entropies"). This variational problem can arise 
in various domains of physics (or biology, economy,...) for different reasons. In any case, it is 
relevant to develop a thermodynamical analogy and use a vocabulary borrowed from thermo- 
dynamics (entropy, temperature, chemical potential, caloric curve, free energy, microcanonical 
and canonical ensembles,...) even if the initial problem giving rise to this variational problem 
is not directly connected to thermodynamics. Thus, we can directly transpose the methods 
developed in the context of ordinary thermodynamics (e.g., Legendre transforms, turning point 
arguments, bifurcations,...) to a new context. For example, in the present study, the max- 
imization of Tsallis entropy at fixed mass and energy is a condition of nonlinear dynamical 
stability for stellar polytropes via the Vlasov equation and for polytropic vortices via the Euler 
equation. On the other hand, the minimization of Tsallis free energy at fixed mass is connected 
to the nonlinear dynamical stability of polytropic stars via the Euler- Jeans equations. It is also 
a condition of thermodynamical stability (in a generalized sense) for self-gravitating Langevin 
particles experiencing anomalous diffusion and a condition of dynamical stability for bacterial 
populations. Although the formalism is the same for all these systems, the results have a very 
different physical interpretation. Our results may also have unexpected applications in other 
domains of physics that we are not aware of. 

On a technical point of view, we have provided the complete equilibrium phase diagram 
of self-gravitating polytropic spheres for arbitrary value of the polytropic index n and space 
dimension D. Our study, generalizing the classical study of Emden jlH] and Chandrasekhar 
[IS], shows how the phase portraits previously reported in the literature (for particular di- 
mensions and particular polytropic indices) connect to each other in the full parameter space. 
From the geometrical structure of the generalized caloric curves, we can immediately deter- 
mine the domains of stability of the polytropic spheres by using the turning point method . 
These stability results have been confirmed by explicitly evaluating the second order varia- 
tions of entropy and free energy. This eigenvalue method provides, in addition, the form of 
the density profile that triggers the instability at the critical points. Interestingly, this study 
can be performed analytically or by using simple graphical constructions in the Milne plane. 
We have found that complete stellar polytropes (with n < n 5 = (D + 2)/(D — 2) if D > 2) 
are stable for D < 4 and unstable for D > 4. On the other hand, complete gaseous poly- 
tropes are stable for D < 2 and for (D > 2, n < n% = D/(D — 2)) and unstable for (D > 2, 
n > n%). Polytropes with index = D/2 correspond to classical white dwarf stars (i.e., a 
self-gravitating Fermi gas at T = 0). They are self-confined only for D < 2(1 + \pX) and they 
are stable only for D < 4. For D > 4, quantum mechanics is not able to prevent gravita- 
tional collapse, even in the non-relativistic regime. In this sense, D = 4 is a critical dimension. 
Therefore, our dimension of space 2<D = 3<4is bounded by two critical dimensions. 
It seems that this remark has never been made before. The description of phase transitions 
in the self-gravitating Fermi gas at non-zero temperature in dimension D will be considered 
in a future article Other possible extensions of our work would be to consider different 
equations of state such as the modified isothermal p = — Tln(l — p/po) associated with an 
"entropy" functional S[p] = — J {plnp + (p — p) ln(po — p)}d D r or the logotropic equation of 
state p = p c [l + A\n(p/p c )} j^I] associated with S[p] = p c A J In pd D r. 
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The concept of generalized thermodynamics is rigorously justified in the case of stochastic 
(Langevin) particles experiencing anomalous diffusion. This happens when the diffusion coef- 
ficient in the Fokker-Planck equation depends on the density of particles while the friction or 
drift is constant. In this paper, we have explicitly studied the nonlinear Smoluchowski-Poisson 
system (for self-gravitating Langevin particles) corresponding to a power-law dependance of 
the diffusion coefficient. This particular situation is connected to Tsallis generalized ther- 
modynamics but more general Fokker-Planck equations can be constructed and studied [20 . 
The connection between thermodynamical and dynamical stability for this type of general- 
ized Fokker-Planck equations has been established in the general case in [20] ■ The nonlinear 
Smoluchowski-Poisson system can have physical applications for the chemotaxis of bacterial 
populations. The collapse and aggregation of bacterial populations is similar, in some respect, 
to the phenomenon of core collapse in globular clusters (or to the Jeans instability in molecular 
clouds) and the neglect of inertia is rigorously justified in biology at variance with astrophysics. 
In addition, biological systems are likely to experience anomalous diffusion so that the NSP 
system can provide an interesting and relevant model for the problem of chemotaxis. We have 
shown that the solutions of the NSP system can either converge toward a complete polytrope 
or an incomplete polytrope restrained by the box, or lead to a situation of collapse. The deter- 
mination of the scaling exponent a in the microcanonical ensemble (constant energy) is difficult 
due to the extremely slow entry of the system in the scaling regime. However, it seems to be 
given asymptotically by a n = 2n/(n — 1) (a = 2 for isothermal spheres) as in the canonical 
ensemble (constant temperature). We expect that an exponent a > a n will be selected when 
the temperature is allowed to vary in space and time. This problem will be considered in a 
future study [3CL 
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A Gravitational force in D dimensions 



The gravitational field produced in r by a distribution of N particles with mass m in a space 
of dimension D is 



(203) F = - Gm 



N 

r - r, 



i=l 



r — v, 



\D ' 



For D — 1, the gravitational field created by a particle is independent on the distance. Thus, 
an object located in x experiences a force (by unit of mass) F = Gm(N + — N~), where N + is 
the number of particles in its right (xj > x) and N~ the number of particles in its left (x« < x) . 

The external gravitational field created by a spherically symmetric distribution of matter 
with mass M is 

(204) F = -V$ = -^e r . 

For D 7^ 2, the gravitational potential is 

GM 



(205) $ = 



(D — 2)r D ~ 2 ' 



where the constant of integration has been taken equal to zero (this implies $ = at infinity 
for D > 2). For D = 2, we have, 

(206) $ = GM\n(r/R), 
where we have taken $(/?) = 0. 

B Virial theorem in D dimensions 

We define the Virial of the gravitational force in dimension D by 

(207) V D = J p r • V$ d D r. 

For a spherically symmetric system, the Gauss theorem can be written 

, nnn , d$ GM(r) , /"" 'Dili 

(208) — = M(r) = ^ P S D r D ^dr>. 

Therefore, the Virial is equivalent to 

nR dMGM(r) G [ R dM 2 1 



In D = 2, one has directly 

GM 2 



dr. 



(210) V 2 = 2 

If now D 7^ 2, we obtain after an integration by parts 

, x GM 2 1 . GM(r) 2 , 
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or, using Eq. (j2U8|) . 

Now, using the Poisson equation (j3J), the potential energy can be written 

(213) W = - I p$d D r = — J— / $A$d D r. 

2 7 2SdG J 

Integrating by parts, we obtain 

(214) W = -1-U(R)^(R)S D R D - 1 - [ R ( ^-Yr D -'S D dr 



2SoG [ dr Jo \ dr 

The gravitational force and the gravitational potential at the edge of the box are given by Eqs. 
f!204|) and ()205|) . Introducing these results in Eq. ()214|) and comparing with Eq. ()212|) . we 
obtain 

(215) V D = -(£>- 2) W, D^2. 

By using the Virial tensor method introduced by Chandrasekhar we can show that the 
foregoing relation remains valid if the system is not spherically symmetric. 
If now the system is in hydrostatic equilibrium, we have 

(216) Vp = -pV$. 

Inserting this relation in the Virial ()2()7)1 and integrating by parts, we get 

(217) V D = - (p pr-dS + 2K, 



where we have used K = (D/2) J pd D r. This is the expression of the Virial theorem in its 
general form. Assuming now that p b is uniform on the domain boundary (which is true at least 
for a spherically symmetric system), we have 

(218) jpY-dS = p b jr-dS = p b J V • r d D r = p b DV D R D . 

Therefore, for a spherically symmetric system, the Virial theorem reads 

(219) 2K-V D = p(R)DV D R D , 

where Vd is given by Eqs. (|21(J|) and (|215|) . 

It is interesting to consider a direct application of these results. In D = 2, the Virial theorem 
reads 

CM 2 

(220) 2K - = 2ttR 2 P {R). 
For an isothermal gas K = MT so that 

GM 2 

(221) 2MT — = 2nR 2 p(R). 

From this relation, we conclude that there exists an equilibrium solution with p(R) = at 
T c = GM/4. We therefore recover the critical temperature of an isothermal self-gravitating 
gas in two-dimensions (see, e.g., 0). At T = T c , the density profile is a Dirac peak so that 
p(R) = 0. For T < T c , there is no static solution and the system undergoes a gravitational 
collapse. This collapse was studied in \3. with the Smoluchowski-Poisson system. 
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C Special properties of n 3 / 2 polytropes 



In this Appendix, we consider polytropes with index 113/2 = D/2 in a space of dimension D > 4. 
They correspond to D-dimensional "white dwarf" stars. According to Eq. ()113|h the curve 77(a) 
is extremum for 

, ^ D(D-A) 
(222) u = K D _ 2 - 

It is easy to check that this particular value is also solution of Eq. (jll7)) determining maxima 
of A(a). We conclude therefore that the functions 77(a) and A(a) achieve extremal values for 
the same values of a in the series of equilibria. This implies that the generalized caloric curve 
77(A) of 77,3/2 polytropes displays angular points (see Figs. Hoi and ITSj) . 



D Second order variations of generalized entropy and 
free energy 

According to Eq. (|35|). the variations of entropy up to second order are 



( ),S' = D - 2n { (3 I 5pd D r + 6(3 I pd D r + / 5(35pd D r 



(223) 

On the other hand, according to the polytropic equation of state (|3~U|). we have 
(224) 

From Eqs. (JSJ) and (ETTl) . 

(225) K ~ (3 1 



bp ip {bp) 2 5K SK5p 
p n 2p 2 K K p 



so that to second order 



(226) 

Inserting Eqs. 



5K _D-2n5(3 (D — 2n) (D — An) ( 5(3 
~ 8^ [jT 



K 2n (3 
and (HP) in Eq. 



5S 



D — 2n 



P 



we get 
Pi f P 



ft I ~pd D r+^ j f^P) 2 d D r+^ / P 5(3d D r 



(227) 



2n J p 

Now, the conservation of energy ()34j) implies that 

D 
~2 

in Eq. 



577/ 







pd r 



(228) 



Inserting Eqs. 

D — 2n 



= 5E 



and 



5pd D r + 5p5$d D r + / <$>5pd D v. 



we obtain 



D-2n P J p P + 2n{D-2n) 



i / i,-L-d D r+^- i S l,„l D r 



D(D - An) (5(3f 
8n~ 2 /T 



pd D r + jRgp [ 6 _R pd D r 



2n 



(229) 



+^(3 j 5p5$d D r + (3 j $5pcTr = 0. 



P 
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Subtracting this relation from Eq. ()227|) . we get 



5S = -p<yn / ^pd D r - ~/3 7 / p^-d D r + / pd D v 

J P 2 J p 2 8n p J 

(230) ~p J 5p5§d D r-(3 J §5pd D r. 

Now, to first order, Eq. (|229jl yields 

5(3 An J SSpdPr + f J J 6 -fpd D v 



(231) 



P ' D(D - 2n) J pd D i 



Substituting this relation in Eq. (|23()|) . we find that the second order variations of entropy are 
given by Eq. ()137|) . To compute the second order variations of free energy F = E — TS, we 
can use Eqs. (E27I) and (J22HJ) with 5(3 = 0. This yields Eq. 



E Some useful identities 

In this Appendix, we establish the identities ([147)1 - (|149j) that are needed in the stability analysis 
of Sec. 13.71 Using an integration by parts, we have 



pa pa j / nn+1 \ i 

(232) / e>ee n di = t°±(l—)dt = -±-cPff*\_ 

Jo Jo d£\n + lj n + 1 n + 1 

Using the Lane-Emden equation and integrating by parts, we obtain 

(233) P e n+1 i D ^di = -a D - l 6(a)6\a) + f° V) 2 ^ -1 ^- 
Jo Jo 

Using the relation 



which results from a simple integration by parts, we get 

D-l 
2 

d£ 



(235) / Z D - 1 (6') 2 dt = -a D 9'(a) 2 + 2 j° ^±.(^0' \&d£, 



or, equivalent ly, 



(236) D I £ D - l {e') 2 d£ = a D e'{a) 2 -2 / ^ D e"e'd£. 

Jo Jo 
Using the Lane-Emden equation (|45|). we find that 

(237) (D - 2) f e D_1 {O'fd^ = -a D 9\a) 2 - 2 f £ D 9'B n d£. 

Jo Jo 



We have three equations (|232|) . (|233|) and (|237|) for three unknown integrals. Solving this 
system of algebraic equations and introducing the Milne variables ([70)1 . we obtain the identities 
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F Dynamical stability of gaseous spheres 



In this section, we assume D > 2. According to Eq. (|57|). the energy of a polytropic star at 
equilibrium (u = 0) can be written 

2n 

(238) W = —K + W, 

where K is the kinetic energy and W the potential energy. Now, for a complete polytrope 
{pb = 0), the Virial theorem reads 

(239) 2K + (D -2)W = 0. 
Combining the foregoing relations, we get 

(240) W=\\-^W. 

According to Poincare's theorem, a gaseous star with W > is unstable JH]- For polytropic 
stars, this condition is equivalent to n > n 3 . 

More generally, the internal energy of a mass dm of gas at temperature T is dU = C v dmT. 
Its kinetic energy is dK = ^dmRT = y(C p — C v )dmT where R is the constant of perfect gases 
and C v , C p are the specific heats at constant volume and constant pressure, respectively. Thus, 
we get 

(241 » u= Dfhr) K 

where 7 = C p /C v . For a monoatomic gas, 7 = (D + 2)/D and U = K. Using the Virial theorem 
(J239j) . the total energy of the star W = U + W can be written 

The star is unstable for 7 < 7 cr jt = 2{D — 1)/D. For D = 3, we recover the well-known result 
Icrit = 4/3. For a polytropic gas, we recover the result 7 cr j t = 1 + l/n 3 . 
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